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PREFACE 


The study described herein was authorized by the Officer In Charge of 
Construction (OICC), TRIDENT, Naval Facilities Engineering Command, Department 
of the Navy, St. Marys, Georgia. All elements of the investigation were con- 
ducted at the US Army Engineer Waterways Experiment Station (CEWES) from 
October 1983 to December 1985. The Coastal Engineering Research Center (CERC) 
conducted this study as part of a larger CEWES modeling effort for Kings Bay 
coordinated by the Hydraulics Laboratory (HL), CEWES, with Messrs. William H. 
McAnally, Jr., and Mitch A. Granat, HL, serving as overall CEWES Project 
Managers. Contract monitoring for the study was provided by Messrs. George 
Carpenter, John Randall, and Brian Smith, OICC. 

This report was prepared by Dr. S. Rao Vemulakonda, Project Manager, 
CERC. The study was performed by Drs. Vemulakonda and Norman W. Scheffner, 
Mr. Jeffrey A. Earickson, and Mrs. Lucia W. Chou, Coastal Processes Branch 
(CR-P), CERC. Work was done under direct supervision of Mr. H. Lee Butler, 
Chief, Research Division, and under general supervision of Dr. James R. 
Houston and Mr. Charles C. Calhoun, Jr., Chief and Assistant Chief, CERC, 
respectively. 

The advice of Dr. Houston and Mr. Butler and the assistance of 
Mr. Bruce A. Ebersole, CR-P, and Ms. Mary A. Cialone (formerly of CERC) are 
acknowledged. Hindcast wave information was provided by Dr. Robert E. Jensen, 
Coastal Oceanography Branch, CERC. This report was edited by Ms. Shirley A. J. 
Hanshaw, Information Products Division, Information Technology Laboratory, 
CEWES. 

During report publication COL Dwayne G. Lee, CE, was Commander and 


Director of CEWES. Dr. Robert W. Whalin was Technical Director. 


EXECUTIVE SUMMARY 


St. Marys Inlet is a large jettied tidal inlet through the barrier 
island system of Georgia and Florida. It is located approximately 30 miles 
north of Jacksonville, Florida. The inlet is the main ocean entrance to the 
US Navy Submarine Base at Kings Bay, Georgia. As a part of upgrading the base 
to accommodate Trident submarines, it became necessary to improve the base 
facilities and modify the navigation channels to the interior and exterior of 
the inlet. 

In 1983 the Officer In Charge of Construction (OICC), Trident, requested 
the US Army Engineer Waterways Experiment Station (CEWES) to determine, by 
modeling, the impact of these modifications on hydrodynamics and sedimenta- 
tion. As a result, two studies were undertaken simultaneously. The first, 
performed by the Hydraulics Laboratory of CEWES, used a hybrid modeling ap- 
proach to determine the impact of the interior modifications. A report by 
Granat, et al. (in preparation) describes the results of the investigation. 
The second study was a numerical modeling effort by the Coastal Engineering 
Research Center (CERC) of CEWES to determine the effect of modifications of 
the exterior channels on coastal processes near the inlet, especially channel 
shoaling rates. The report herein describes details of this second study. 

To accomplish the objectives of the second study, CERC employed a system 
of numerical models called Coastal and Inlet Processes (CIP) Numerical Model- 
ing System. The system included four separate numerical models for tides, 
waves, wave-induced currents, and noncohesive sediment (sand) transport. The 
system together with two computational grids--one for existing (base) condi- 
tions and the other for plan conditions--was called Model B in contrast to 
Model A, the hybrid model for the interior. 

To substantiate the validity of the modeling approach and to improve 
accuracy of predicted results, Model B was first verified with available field 
data on tides and sediment transport. The tidal model was verified by using 
field data on tidal elevations and currents taken on 10 November 1982. Good 
verification was obtained by matching model results with observed tidal cur- 
rents over one tidal cycle. The wave climate for an average year in 60-ft 
depth mean low water (mlw) at the project site was obtained from the CEWES 
Wave Information Study based on a 20-year hindcast. This information was dis- 


cretized into 79 different monochromatic waves. These waves were propagated 


to the shore using the wave model and wave conditions were determined every- 
where over the study area. The wave-induced current model used the wave in- 
formation to determine wave-induced currents over the study area. The sedi- 
ment transport model used the results of the other three models to determine 
sediment transport. It was verified by comparing its predictions on naviga- 
tion channel shoaling rates with shoaling rates computed from channel surveys 
taken during the period November 1980 to December 1981. There was good agree- 
ment with respect to both trends and magnitudes. 

Model B was next used to determine base conditions corresponding to 
trapezoidal entrance and offshore channels with a bottom width of 400 ft, a 
project depth of 40 ft mlw and side slopes of 4H:1V. The sediment transport 
model results were obtained in terms of channel shoaling rates (ft/year) along 
the channel. The results were similar to those obtained during verification. 
In both cases, there was deposition outside the jetty tips. It changed to 
erosion inside the jetty tips because of circulation due to wave-induced cur- 
rents. The heaviest deposition rates were predicted near the jetty tips. 

This is the area where the channel cuts through the offshore bar and where 
serious shoaling problems were encountered in the field for base conditions. 
On the basis of the numerical results, the yearly channel shoaling volume 
between sta -80+00 and sta 325+00 was predicted to be 475,000 cu yd/year. 
This value was comparable to within +25 percent of the yearly maintenance 
dredging volumes recorded by the US Army Engineer District, Jacksonville. 

Model B tested only one plan condition which was called Plan 1. The 
plan was to (a) widen the navigation channel to 500 ft, with the widening tak- 
ing place on the north side of the present entrance and offshore channels; 

(b) extend the channel on the ocean side, with the extension being at an angle 
of 20 deg south of the present channel center line at sta -97+76 approxi- 
mately; and (c) deepen the channel to -49 ft mlw (46-ft project depth plus 
3-ft advance maintenance). The channel was to have a trapezoidal cross sec- 
tion with side slopes of 3H:1V. At the request of OICC, TRIDENT, it was 
assumed during testing that the landward 1,000 ft of the south jetty would be 
made sand tight simultaneously. 

In view of the urgent needs expressed by OICC for Plan 1 results from 
Model B for design of entrance and offshore channels, the wave and wave- 
induced current models were not rerun for Plan 1 as originally planned. Only 


the tidal and sediment transport models were rerun. The results of the tide 


model showed the effects of Plan 1 to be mainly local and caused by sand- 
tightening of the south jetty. Tidal velocities at the end of the jetties and 
at Tide Gage 1 near the south jetty increased by approximately 10 percent. 
There were no significant changes in velocities at the throat of the inlet. 
The extension of the navigation channel at the seaward end produced almost no 
effect upon the tidal current patterns. 

The results of the sediment transport model indicated an increase in 
both deposition rates and erosion rates from base to Plan 1 for the channel 
reach between sta -97+76 to sta 325100. Model results for the reach between 
sta 325+00 and sta 399+74 for base and Plan 1 were suspect since quantitative 
field information on sedimentation was not available for this reach, the 
bathymetry used in Model B was not the latest, and the grain size of the sedi- 
ment observed in this reach was much larger than that of the sediment every- 
where else in the study area. For the channel reach between sta —80+00 and 
sta 325+00, the shoaling volume for Plan 1 was estimated from Model B results 
to be approximately 788,000 cu yd/year or an increase of approximately 66 per- 
cent from base. Finally, based on Model B results for Plan 1 and all other 
available information, recommendations on advance maintenance dredging were 
made for various reaches of the channel for use in channel design. Based on 
modeling limitations, the accuracy of Model B sediment transport results is 


estimated to be within +25 percent. 
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CONVERSION FACTORS, NON-SI TO SI (METRIC) 
UNITS OF MEASUREMENT 


Non-SI units of measurement used in this report can be converted to SI 


(metric) units as follows: 


Multiply By To Obtain 

cubic feet per second 0.0929 cubic metres per second 

per foot per metre 
cubic yards per year 0.7646 cubic metres per year 
feet 0.3048 metres 
inches 0.0254 metres 
miles (US nautical) 1,852 metres 
miles (US statute) 1,609 metres 
pounds (force) per foot 14.5932 newtons per metre 
pounds per second 0.4536 kilograms per second 
pounds (force) per 47.88 pascals 


square foot 


slugs per cubic foot 515.4 kilograms per cubic metre 


KINGS BAY COASTAL PROCESSES NUMERICAL MODEL 


PART I: INTRODUCTION 


Background 


1. St. Marys Inlet is a large jettied tidal inlet through the barrier 
island system of Georgia and Florida. It is the main entrance to Kings Bay 
Naval Submarine Base located at Kings Bay, Georgia. The inlet is located 
approximately 30 miles* north of Jacksonville, Florida (Figure 1). The 
Georgia-Florida state line runs through the inlet. To the north of the inlet 
is Cumberland Island administered by the National Park Service, and to the 
south is Amelia Island. Fort Clinch State Park surrounding historic Fort 
Clinch is located on Amelia Island. 

2. At present, Kings Bay is home to Poseidon-class submarines. The 
present entrance and offshore channels are trapezoidal in cross section with a 
bottom width of 400 ft, a project depth of 40 ft mean low water (mlw**), and 
side slopes of 4H:1V. As a part of the upgrading of the submarine base to 
receive the larger Trident-class submarines, it became necessary to widen and 
deepen both the interior and exterior navigation channels. Simultaneously, it 
is proposed to sand-tighten a 1,500-ft segment of the south jetty. This study 
is mainly concerned with the exterior (entrance and offshore) channels. Here- 
after the term "entrance channel" will be used to refer to the part of the 
exterior navigation channel between the jetties, and the term "offshore chan- 
nel" will be used to denote the part of the navigation channel offshore of 
jetty tips. A companion study (Granat, et al., in preparation) considers the 


interior channels. 


Purpose 


3. The purpose of this study is to determine the effect of the project 


on coastal processes near St. Marys Inlet. The processes studied include 


* A table of factors for converting non-SI to SI (metric) units is presented 
on page 6. 
**k Abbreviations and acronyms are listed in Appendix B. 
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tides, waves, wave-induced currents, and sediment transport. Of special 
interest is the determination of shoaling rates in the navigation channel for 
existing (base) and plan conditions. Based on these shoaling rates, recom- 
mendations will be made to the project channel designers on the required 
yearly advance maintenance dredging for various reaches of the channel. To 
accomplish the objectives of the study, a system of numerical models called 
Coastal and Inlet Processes (CIP) Numerical Modeling System, which includes 
separate numerical models for the coastal processes mentioned above, will be 
employed. The system will take into account the effect of the two jetties on 
St. Marys Inlet. Models of the system will be calibrated and verified with 
available field data as far as possible and used to study existing conditions 
as well as planned conditions to determine the effect of the project on 


coastal processes. 


PART II: THE CIP NUMERICAL MODELING SYSTEM 


Numerical Modeling 


4. Before details of the CIP modeling system are presented, a brief in- 
troduction to certain aspects of numerical models is in order. Generally, the 
physical variables of practical interest such as surface elevation, velocity, 
wave height, and sediment transport rate vary continuously in space and time. 
On the basis of the physics of the particular process, the variables are de- 
scribed by differential equations. In numerical modeling, the differential 
equations are replaced by difference equations involving finite differences in 
space and time. Thus, a numerical model considers values of the variables at 
discrete points in space and time and solves for the values of the variables 
by numerical techniques. 

5. Numerical models are classified on the basis of variation of the 
dependent variables in space and time. If the dependent variable is a func-— 
tion only of one coordinate, then we have a one-dimensional model. For exam— 
ple, the average velocity in a river cross section may be a function only of 
distance along the river, and we can describe the flow using a one-dimensional 
model. If the dependent variable is a function of two coordinates, then we 
have a two-dimensional model. For instance, tidal elevations and currents in 
a shallow bay may be a function only of the two horizontal coordinates, and 
the tidal hydrodynamics can be described by a vertically averaged two- 
dimensional model. If the dependent variable is independent of time, a steady 
model is applicable; whereas if the variable varies with time, an unsteady 
model is needed. 

6. It should be recognized that numerical models are only as good as 
the physics that goes into them and are in general approximations to physical 
reality. In recent years, numerical models have become standard tools to 
answer questions connected with engineering projects and have replaced tradi- 
tional physical hydraulic models for studies involving tidal hydraulics, wave 
transformation, etc. They are the only feasible tools available for analyzing 
certain phenomena such as sediment transport under the combined action of 
tides and waves, wind-generated flows, etc. They have the following advan- 
tage. Once a numerical model has been calibrated and verified for a given 


project area for a given set of conditions, it can predict, within a 


10 


reasonable degree of reliability, the physical processes under a different set 
of conditions, provided the latter set is not radically different from the 
first. Thus the model is usually calibrated and verified for previous or 
existing conditions and used to predict future plan conditions. 

7. In the study described herein, the coastal processes for St. Marys 
Inlet and the surrounding area of the Atlantic Ocean were modeled using the 
CIP numerical modeling system on two computational grids. The system includes 
the US Army Engineer Waterways Experiment Station (CEWES) Implicit Flooding 
Model (WIFM) for tides, the Regional Coastal Processes Wave Propagation Model 
(RCPWAVE) for waves, the model CURRENT for wave-induced currents, and a sedi- 
ment transport model for transport of noncohesive sediments due to the com- 
bined action of tides, waves, and wave-induced currents. All four models 
generally used the same computational grid for a given set of conditions (base 
or plan). The following paragraphs highlight the important features of the 
computational grids and the four computer models used in this study. For con- 
venience the numerical modeling system, together with the computational grids, 
was referred to as Model B in contrast to Model A, a hybrid model used for 


studying the region interior to the inlet. 


Computational Grids 


8. The models described in this report use the finite difference method 
for computations. In order to cover a large region but still maintain high 
resolution in desired areas, the models use a smoothly varying grid that 
allows cells to be small in certain areas (e.g., surf zone or inlet) and large 
in others (e.g., ocean or sound). A piecewise reversible transformation 
(analogous to that used by Wanstrath (1977)) is used independently in the x- 
and y-directions to map the variable grid into a uniform grid used in the com- 


putational space (Figure 2). The transformation has the following form: 


ce 

x=a_+ boo" (1) 
c 

y=a_+bD an (2) 


11 


Wns & 5 WD 56 © 5» & 5 WD 5 el Cio Are Alelslinrehay COMSIERMIES iter 
P P Pp q q q 


regions p and q in the x- and y-directions, respectively, and Oy and a 


are coordinates in the computational space. This transformation allows all 


2 


derivatives to be centered in the computational space. Many stability prob- 
lems commonly occurring in variable grid schemes are eliminated when using 
this transformation since the grid in real space varies smoothly and the co- 
ordinates and their first derivatives are continuous at the boundaries between 


regions. 
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Figure 2. An example of variable grid 


9. The partial differential equations governing the different processes 
are solved by finite difference integration on a grid of spatial points. A 
right-handed coordinate system is used with the x-coordinate increasing in the 
offshore direction and the y-coordinate increasing along the shoreline with 
the ocean to the right. The partial derivative of an arbitrary variable s 


in region p is 


OS I os 
aan e oas (3) 


where 


(4) 


1= 
i] 
ion 
Q 
eg 


* For convenience, symbols and abbreviations are listed in the Notation 
(Appendix A). 
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Similarly 


as _ i ds 
Smet 3a, (5) 


where 


ec -l 


eet otras q 
uy 205 Dean (6) 


Q 


If the grid in the x-, y-coordinate system is to have constant grid spacing, 
all values of Hy. and ne will be constant (1 if Aa, = Ax and Aa, = Ay). 
The constants a. He te |b eat eee Cedar iment: emp’? aniaes| ome area 8 01S e for all the regions and 


the values of uy pean Me a a4 Ne ene and centers are determined 
using an interactive computer program called MAPIT. A plotting program called 
CMSGRID is used to plot the grids to a suitable scale for overlaying nautical 
charts such as those of the National Ocean Service (NOS). 

10. Figures 3 and 4 show the two computational grids used to model 
St. Marys Inlet and surrounding areas. Both grids extend for 141,670 ft 
(23.3 nautical miles), in 90 grid cells, east-west and for 60,000 ft (9.9 nau- 
tical miles), in 50 cells, north-south. The grids are oriented so that one 
row of cells aligns with the navigation channel. The difference between the 
two grids lies in the mapping of 17 rows of cells in the area which covers 
St. Marys Inlet. Grid 1 (Figure 3) models the base condition of a 400-ft-wide 
navigation channel. The smallest cell, located near Ft. Clinch, has dimen- 
sions of 275 ft by 310 ft; while the largest cell, located at the eastern 
boundary of the grid, is 5,200 ft by 4,170 ft. The 25th cell row from the 
north edge of the grid overlays the channel. Grid 2 (Figure 4) models the 
plan condition of a 500-ft-wide channel. The smallest cell in Grid 2, also 
near Ft. Clinch, is 229 ft by 310 ft, while the largest cell size remains the 
same. The 17 rows of cells in the area of St. Marys Inlet have been remapped 
for Grid 2 so the 25th row again overlays the channel. This adjustment is 
accomplished by mapping slightly larger cells just north of the channel and 
slightly smaller cells just south of the channel. The mapping differences 
between Grids 1 and 2 are so minor that no bathymetic changes are required 


other than those in areas affected by dredging for the plan conditions. 
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11. Due to differences in required boundary conditions between the 
various models, and for reasons of economy, the first 17 columns of cells 
along the west side of the grids (the interior area) were not used in RCPWAVE, 
CURRENT, and the sediment transport model. This area is covered by Model A. 
WIFM required the additional cells in the interior area in order to use the 
prototype data available there for boundary conditions and to accurately simu- 


late the complex tidal currents in the inlet. 


The Tidal Simulation Model, WIFM 


12. WIFM is a general long wave model which can be used for simulation 
of tides, storm surges, tsunamis, etc. It allows flooding and drying of land 
cells near the shoreline. It is a depth-averaged model so that variations in 
the vertical direction are averaged in the model. It is used in the present 
study to determine tidal elevations and velocities in the two horizontal 
coordinate directions. The following description of WIFM is extracted from a 
report by Leenknecht, Earickson, and Butler (1984). 

Equations of motion 
13. The hydrodynamics of the numerical model WIFM are derived from the 


Navier-Stokes equations in a Cartesian coordinate system (Figure 5). The long 


WATER SURFACE 


BENCHMARK CATUM 


Figure 5. Coordinate system for WIFM 
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wave approximations of small vertical accelerations and a homogenous fluid 


yield the following vertically integrated (depth-averaged) two-dimensional 


equations of continuity and momentum: 


Continuity 
an 0 0 
ee ee (weal) ay (wel) = (7) 
Momentum 
YP 
du du du a gu 2 
= —= —- — - + + 
eee ee ge Oa) tala ea VD) 
Cid 
Zz 
2 2 
= 6 gu,au cist neh e—iai(Q) (8) 
2 2 x 
ox oy 
1/2 
ov OV OV Q) gv 2 
— —+y + +e — = a + 
are iP Wl rev By fu + g 9 (n n,) Pa (u va) 
Cd 
Zz 
2 2 
cae over oy ER Ee (9) 
ox oy y 
where 
nN = water surface elevation above datum 
t = time 
u,v = velocities in the x- and y-directions 
d =n +h, the total water depth 
h = local still-water depth 
R = rate of water volume change in the system due to rainfall or 
evaporation 
= Coriolis parameter 
g = acceleration due to gravity 
Cc = Chezy coefficient for bottom friction 
€ = eddy viscosity coefficient 
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The variable is accounts for hydrostatic water elevations due to atmospheric 
pressure differences, and Es and By represent external forces such as wind 
stress. 
Numerical method 

14. The alternating-direction-implicit (ADI) method has been used by 
Leendertse (1970) and others to solve the two-dimensional equations of motion. 
When the advective terms are included in the momentum equations (Equa- 
tions 8, 9) the ADI method has encountered stability problems. Weare (1976) 
indicates that the problems arise from approximating advective terms with 
one-sided differences in time and suggests the use of a centered scheme with 
three time-levels. WIFM employs a centered stabilizing-correction (SC) scheme 
which is second-order accurate in space and time, and boundary conditions can 
be formulated to the same order of accuracy. A brief development of the SC 
scheme is presented in the following paragraphs. Note that n and h are 
defined at the cell center and u and v _ at the cell faces. 

15. The linearized equations of motion can be written in matrix form 


as: 
Wig AUL BUY = 0 (10) 


where 


The SC scheme for solving Equation 10 is 
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k-1 
k= = oO 
(1 + ) U (1 ne 22.) U Ci) 
@eaqj) tf emer ee (12) 
y y 
where 
1 At 
i iF 2 Ax AS 
1 At 
A == — BS 
SZ BY Oy 


The quantities Oy. and 8 are centered difference operators, and the super- 
script k indicates time-level. The starred quantities can be considered 
intermediate values between the k and k+l time-levels. 

16. The first step in the SC procedure computationally sweeps the grid 
in the x-direction, with the second step sweeping in the y-direction. Com- 
pleting both sweeps constitutes a full time-step, advancing the solution from 
the k time-levell’ to the kt! time-level- The form of the difference equa— 


tions for the x-sweep is given by 


a qe eo we") + xs, 6. GER a) & = S (ta) = 0 (13) 
— (ey = oy 2 a= CG = qo) 30 (14) 
— Cle eye ts oF GD) = @ (15) 
and the y-sweep by 
Go eae as 7 ase (16) 
ysth oe (17) 
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1 k+l g k+1 k-1 
— = * = = 
(v ve) + DAy o, (n n ) 0 (18) 


17. Noting that v* in Equation 15 is only a function of previously 
computed variables at the k-l time-level, its substitution into Equation 18 
and the substitution of u* (Equation 17) into Equations 13 and 14 yield the 


simplified forms 


x—Sweep 
se (ne TH) te a a + uh Tay + is s, Go) sO ~ C9) 
eee aha Lal lt on Meneame ee 

y-sweep 
Reena as ae a pace acc mle i ae 
eo ie thn i 8, sh gia, 2565 


18. The details of applying the SC scheme to Equations 7-9 can be found 
in a report by Butler (in preparation). The diffusion terms of Equations 8 
and 9 are also represented with time-centered approximations. The inclusion 
of diffusion terms contributes to the numerical stability of the scheme 
(Vreugdenhil 1973) and serves to somewhat account for turbulent momentum 
dissipation at the larger scales. While the resulting finite difference forms 
of Equations 7-9 appear cumbersome, they are efficient to solve. Application 
of the appropriate equation to one row or column of the grid (the "sweeping" 
process) results in a system of linear algebraic equations whose coefficient 
matrix is tridiagonal. Tridiagonal matrix problems can be solved directly, 
without the cost and effort of matrix inversion. 

19. Apart from Courant number considerations, the computational time- 
step for the SC scheme in WIFM is largely governed by simple mass and momentum 


conservation principles. The maximum time-step for a problem is characterized 


by 
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At = v (23) 


where V is the largest flow velocity to be encountered at a cell with its 
smallest side length AS . The parameter n is of order 1. Therefore, the 
time-step is constrained by the smallest cell width which contains the highest 
flow velocity. In physical terms, Equation 23 requires that the flow cannot 
move substantially farther than one cell width in one time-step. 

Boundary conditions 

20. WIFM allows a variety of boundary conditions to be specified, which 
can be classified into three groups: open boundaries, land-water boundaries, 
and thin-wall barriers. 

21. Open boundaries. When the edge of the computational grid is 
defined as water, such as a seaward boundary or a channel exiting the grid, 
either the water elevation or the flow velocities can be specified as an open 
boundary-condition. This information can be input to WIFM as tabular data, or 
constituent tides can be calculated within the model during the time-stepping 
process. 

22. Land-water boundaries. WIFM allows land-water boundaries to be 
either fixed or variable to account for flooding in low-lying terrain. Fixed 
boundaries specify a no-flow condition at the cell face between land and 
water. The position of a variable boundary is determined by the relationship 
of the water elevation at a "wet" cell to the land elevation at a neighboring 
"dry" cell. Once a water elevation rises above the level of adjacent land 
height, water is initially moved onto the "dry" cell by using a broad-crested 
weir formula (Reid and Bodine 1968). When the water level on the dry cell 
exceeds some small value, the boundary face is treated as open, and computa- 
tions for n ,u , and v are made at the now "wet" cell. Drying is the 
inverse process, and mass is conserved in these procedures. 

23. Thin-wall barriers. These barriers are defined along cell faces 
and are of three types: exposed, submerged, and overtopping. Exposed 
barriers allow no flow across a cell face. Submerged barriers control flow 
across a cell face by using a time-dependent friction coefficient. Overtop- 
ping barriers are dynamic. They can be completely exposed, completely sub- 
merged, or they can act as broad-crested weirs. The barrier character is 


determined by its height and the water elevations in the two adjoining cells. 
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The Wave Model, RCPWAVE 


24. The RCPWAVE model is a linear short-wave model which considers 
transformation of surface gravity waves in shallow water including the pro- 
cesses of shoaling, refraction, and diffraction due to bathymetry and allows 
for wave breaking and decay within the surf zone (the region shoreward of the 
breaker line). Unlike traditional wave-ray tracing methods, the model uses a 
rectilinear grid so that model output in the form of wave height, direction, 
and wave number is available at the centers of the grid cells. This avail- 
ability is highly advantageous since the information can be used directly as 
input to the wave-induced current and sediment transport models, and the prob- 
lem of caustics due to crossing of wave rays is avoided. The description of 
RCPWAVE that follows is extracted from a report by Ebersole, Cialone, and 
Prater (1986). 

25. Berkhoff (1972 and 1976) derived an elliptic equation approximating 
the complete wave transformation process for linear waves over an arbitrary 
bathymetry constrained only to have mild bottom slopes (thus the designation 
mild slope equation (Smith and Sprinks 1975)). The mild slope equation can be 


expressed in the following form: 


Oia, Oa 2s fy SON a 2 Te, o 
x (cc, at ay (:c, oe Gone ¢=0 (24) 


where 


¢(x,y) = complex velocity potential 
27 
Oo = wave angular frequency = rT 
T = wave period 


c(x,y) = wave celerity = 


tam Ke} 


= pe, lee 
eC) group velocity = ak 


k(x,y) = wave number given by the dispersion relation 


of = gk tanh(kh) (25) 


26. Numerical solution of this equation for the velocity potential 


field is an effective means for solving the complete wave propagation problem. 


All 


The equation can be solved using either finite element (for example, Berkhoff 
1972, Houston 1981) or finite difference methods (for example, William, 
Darbyshire, and Holmes 1980). Since transmission and reflection boundary con- 
ditions are easily implemented into these solution schemes, this approach is a 
popular one for modeling tsunami propagation and for solving problems involv- 
ing the response of harbors to short and long waves. This method becomes 
computationally infeasible for large scale, open coast, short-wave problems 
because of its great expense. 

27. The model RCPWAVE is an alternative approach for solving the open 
coast wave propagation problem. It addresses the processes of refraction and 
diffraction and can be applied to a large region quite economically. The 
model also contains an algorithm which estimates wave conditions inside the 
surf zone. This wave breaking model is an extension of the work of Dally, 
Dean, and Dalrymple (1984) to two horizontal dimensions. 

Wave transformation outside 
the surf zone: theoretical basis 
28. The velocity potential function for linear, monochromatic, plane 


waves can be represented by the following expression: 
¢=ae (26) 


where 


& H(x,y) 
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a(x,y) = wave amplitude function equal to 

H(x,y) = wave height 

s(x,y) = wave phase function 
Here the velocity potential function describes only the forward scattered wave 
field. No considerations are given to wave reflections. By substituting this 
expression for the velocity potential into Equation 24 and solving the real 
and imaginary parts separately, two equations can be derived (Berkhoff 
(1976)), namely, 

2 2 
coma’ + -— + ay (v2 ° vee,)) oF a - lnalj = 0 (27) 


ox dy ace 


Ve (a°ce, Vs) = (28) 
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where the symbol V denotes the horizontal gradient operator. 

29. Together, these equations describe the combined refraction and 
diffraction process. Diffraction is often erroneously described as the 
propagation of energy along wave crests which are defined to be perpendicular 
to the wave phase function gradient Vs . Equation 28 shows energy is still 
propagated in a direction perpendicular to the wave crest. Diffractive 
effects do change the phase function as a result of significant gradients and 
curvatures of the wave height. These changes cause the local wave direction 
to vary. If diffractive effects are neglected, Equations 27 and 28 reduce to 
those describing pure refraction in which the wave number represents the mag- 
nitude of the phase function gradient. 

30. Linear wave theory assumes irrotationality of the wave phase 


function gradient. This property can be expressed mathematically as 
Vx (Vs) = 0 (29) 
The phase function gradient can be written in vector notation as 


> 


> 
Vs = |vs | cos 6@it [Vs | sin 6 j (30) 


> > 
where i and j are unit vectors in the x- and y-directions, respectively, 


and 0(x,y) is the local wave direction. Equations 29 and 30 can be combined 


to yield the following expression: 


0) 0) 
a (|vs| sin 6) - ay (|vs| cos 6) = 0 (31) 


If the magnitude of the wave phase function is known, local wave angles can be 
calculated from Equation 31. Similarly, Equation 30 can be substituted into 


Equation 28 to yield 


p) 2 p) 2 . 
ae (a cc, |vs| cos 6) + aH (a cc, |Vs| sin 9) = 0 (32) 


This form of the energy equation can be solved for the wave amplitude function 


a once the wave phase characteristics Vs and 6 are known. The wave 
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height can be determined and is proportional to the amplitude function since 
wave frequency is constant. 

31. Equations 27, 31, and 32 along with the dispersion relation 
describe the combined refraction and diffraction process for linear plane 
waves subject to the restrictions that the bottom slopes are small, wave 
reflections are negligible, and any energy losses are very small and can be 
neglected. These equations are assumed to be valid outside the surf zone. 

The numerical solution scheme used to solve these equations is presented in 
the next section. 

Wave transformation outside 

the surf zone: numerical solution 

32. The three governing equations (27, 31, and 32) are solved using 
numerical methods. Partial derivatives within the equations are approximated 
using finite difference operators. Finite difference solution methods require 
the construction of a computational grid system or mesh. Solution accuracy is 
directly related to resolution within the grid system. Discussions throughout 
this section refer only to grid systems comprised of constant sized, rectangu- 
lar cells. RCPWAVE is capable of computing solutions on variably sized, rec- 
tilinear grid systems. 

33. Figure 6 shows nine rectangular cells which make up a small part of 
a larger mesh. Each cell has a length equal to Ax in the x-direction and 
Ay in the y-direction. The maximum values of i and j are M and N, 
respectively. All variables which vary as a function of space are defined at 
the cell centers (see Ebersole, Cialone, and Prater (1986) for details of the 
finite difference procedure used). 

34. Model input includes values of the deepwater height Hy » direction 
on » and period T of waves to be simulated. It also includes specification 
of the bottom bathymetry throughout the grid. The wave number, which is 
related to the wave period and the local water depth through the dispersion 
relation, is computed at every cell. It is used as an initial guess for the 
magnitude of the wave phase function gradient. The wave celerity c and the 
group velocity e are functions of the wave period, wave number, and water 


depth. Therefore these variables can be calculated at each cell. 
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Figure 6. Definition of coordinate system and grid cell 
conventions used in RCPWAVE 


35. From Snell's law, 


aa waren (88) 


where cy is the deepwater wave celerity (defined to be gy, an estimate of 


the local wave angle is obtained everywhere. This estimate assumes that the 
bottom contours are parallel with the y-axis. If the bottom bathymetric 
contours make a known nonzero angle with the y-axis, a better first guess for 


the wave angles can be made. The new approximation is 


sin(®@ - 6 ) 
Dee Oe Ce a6 (34) 


ZS) 


"The local wave angle, deepwater wave 


where 6, defines the "contour angle.' 
angle, and contour angle follow the angle convention shown in Figure 7. The 
contour angle is an input parameter for RCPWAVE. 


y- AXIS 


POSITIVE 0, 


NEGATIVE 6, 


6 
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@ =LOCAL WAVE ANGLE 
6, =OFFSHORE CONTOUR ANGLE 


Figure 7. Definition of angle conventions used in RCPWAVE 


36. Wave heights at each cell are estimated as the product of the deep- 
water wave height, a shoaling coefficient we and a refraction coefficient 


K  , thus 
r 


H = HK re (35) 
where 
cos on 1/2 
Kr a cos 8 (36) 
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and 


1 1/2 


K = —— 0.0... .C.C.:.:.:. se _  - 
s 2kh 
E + sarc | tanh(kh) 


(37) 


The dispersion relation, Snell's law, and this simple estimator of the wave 
height allow an initial guess to be made for the variables of interest 
throughout the grid system. 

37. The solution scheme implements the following marching procedure 
once initial guesses for the variables of interest have been made. Starting 
at the offshore row eeionaced by i=M-3 , Equations 31 and 32 are used to 
compute wave angles and then heights along the entire row (from j=2 to j=N-1). 
Wave height is used interchangeably with amplitude function since one is 
directly proportional to the other. 

38. Wave angles and heights along a given row are solved for itera- 
tively because of the implicit differencing formulation used. Calculations of 
the wave angle (actually the sine of the wave angle) and the wave amplitude 
function are repeated until the average change (along a row) in each variable 
from one iteration to the next is less than some tolerance. These convergence 
criteria, 0.0005 for sines of the wave angles and 0.001 ft (or a metric equiv- 
alent) for wave heights, are suggested values for prototype applications. 

39. This solution considers only refraction since the wave number k 
is used as an estimate of the magnitude of the phase function gradient. Equa- 
tion 27 is then used to compute the true magnitude of the wave phase gradient. 
This "new wave number" accounts for the effects of diffraction. Backwards 
differences are used to approximate the x-derivatives because they only re- 
quire information which has already been computed. Next, Equations 31 and 32 
are again solved in order to compute the wave angles and heights using these 
new wave numbers. This procedure is repeated along the row under considera- 
tion until the change in new wave number, from one iteration to the next, is 
less than 0.5 percent of the newly computed value. This condition must be met 
at each cell along the row. As a row of new wave numbers is computed, the 
values are filtered in the y-direction using the method of Sheng, Segur, and 
Lewellen (1978). This filter removes cell-to-cell oscillations introduced as 
a result of the differencing scheme used to compute the new wave numbers. 


Row-by-row marching proceeds until solutions are computed along row i=2. 


ZA] 


40. Lateral boundary conditions for a row are specified at the conclu- 
sion of calculations for that row. The value of all variables at cells j=N 
and j=l are set equal to their values at cells j=N-l1 and j=2 , respec- 
tively. This boundary condition implies that the change in the variable in 
the y-direction is zero. The condition is most valid when the bathymetric 
contours are nearly straight and parallel to the y-axis. For this reason the 
grid is oriented so that the y-axis is nearly parallel to bottom contours 
along the lateral boundaries. 

41. Boundary conditions along the offshore boundary of the grid are 
used to initiate the shoreward marching algorithm. They are computed from 
deepwater wave input supplied by the user along with the following assumption. 
Bottom contours extending from the offshore grid row (i=M) out to deep water 
are assumed to be straight and parallel to a line making an angle of 00 with 
the y-axis. In other words, Snell's law is assumed to be valid from deep 
water to the outer boundary of the grid system. No inshore boundary condi- 
tions (along row i=l) are required because of the forward marching solution 
scheme. 

Wave transformation 
inside the surf zone 

42. Waves approaching the very nearshore zone tend to steepen and 
eventually break because of decreasing water depths. Shoreward of this break- 
ing point dissipative energy losses due to turbulence strongly influence the 
wave height. Linear theory does not allow for prediction of the breaker loca- 
tion nor for wave transformation across the surf zone. Instead, empirical and 
approximate methods must be used to describe the breaking process. 

43. The first aspect to consider in surf zone transformation of waves 
is incipient wave breaking. RCPWAVE uses the following criterion of Weggel 
(1972): 


H, = ——— (38) 


where 


Hy = breaking wave height 
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1.56 


fs te 
@ a ee) 
m = bottom slope 
hy = water depth at breaking 
A018. @oor™) 


because it accounts for bottom slope and wave period. 

44. Once the incipient breaking point is defined, a mechanism is needed 
to transform the breaking wave across the surf zone. The transformation 
algorithm selected for use in RCPWAVE (Dally, Dean, and Dalrymple 1984) uses 
an energy flux basis. Through analogy with energy loss in a hydraulic jump in 
a channel, the following equation is postulated for one-dimensional transfor- 


mation of waves advancing in the -x direction: 


acre ss 
Goo i fee, m Gey) | on 
where 
Ec_ = energy flux associated with the breaking wave 
K = rate of energy dissipation coefficient (set equal to 0.2 in 
RCPWAVE) 
(Ec_) = stable level of energy flux that the transformation process 


8 s seeks to attain 


The right side of Equation 39 is simply a dissipation term. The subscript s 
is used to denote the stable level of a variable. Substituting the linear 
wave theory estimate for E (E = 0.125 a) into Equation 39 results in the 


following expression: 


2 
d(H c_) 
Be Ee cu He 
dx hl? sg (x a) (40) 


45. Various field (Thornton and Guza 1982) and laboratory (Horikawa and 
Kuo 1966) experiments have shown that, well into the surf zone, the wave 
height tends toward a stable value which is proportional to the local water 


depth. This relationship can be expressed as 
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H_ = Th (41) 


where 


ia 
I 


stable wave height 
T = proportionality coefficient (set equal to 0.4 in RCPWAVE) 


Equation 40 can now be rewritten as 
K 2 2. 2 

=-|]Hec - Pla & =D 42 
h g ( a oe 


46. This surf zone wave transformation model, extended to two dimen- 
sions, can be incorporated into the conservation of wave energy equation 
(Equation 28) by simply adding a dissipation term D to the right side. The 
function D must now represent dissipation in the direction of wave prop- 
agation. Also for dimensional consistency, the term D must be multiplied by 
the wave celerity and the magnitude of the wave phase gradient, and the wave 
height must be replaced by the wave amplitude function. In vector notation, 


the energy equation becomes 


Vw (ace, Vs) = = ace, |Vs| - (&) rh? ce, |¥5| 5 (43) 
This equation can be thought of as being valid both inside and outside the 
surf zone. Outside, the coefficient k is zero, and the equation reduces to 
Equation 28. 

47. All discussion relating to wave transformation within the surf zone 
up to this point has addressed the problem of determining wave heights. The 
problem of wave phase must be addressed also. Diffraction effects are assumed 
to be negligible inside the surf zone. Therefore, the wave number k is as- 
sumed to accurately represent the magnitude of the wave phase function gradi- 
ent. The linear wave theory assumption that the waves are irrotational also 
will be assumed to remain valid inside the surf zone. Consequently, wave 
angles are computed in the same manner as outside the surf zone. Details con- 
cerning the numerical solution inside the surf zone can be found in Ebersole, 


Cialone, and Prater (1986). 
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The Wave-Induced Current Model, CURRENT 


48. When waves break and decay in the surf zone, in general they induce 
currents in the longshore and cross-shore directions and changes in the mean 
water level. These currents play a major role in the movement of sediment in 
the nearshore. They are computed using the model CURRENT. 

Equations of motion 

49, The hydrodynamic equations used in the model for wave-induced 
currents may be derived from the Navier-Stokes equations (for details, see 
Phillips 1969 and Ebersole 1980). It is assumed in the derivation that the 
fluid is homogeneous and incompressible, and the vertical accelerations are 
negligible so that the pressure distribution is hydrostatic. By vertically 
integrating the three-dimensional form of the equations and applying appropri- 
ate boundary conditions, the depth-averaged two-dimensional form of the equa- 
tions of motion and continuity are obtained. These equations are derived by 
time-averaging over a time interval corresponding to the period of the waves. 


Referring to a Cartesian coordinate scheme (Figure 8), these are: 


Momentum 


= chs) 0s OT 
OW a g§ Magy Beg Mos g “ul x, Sy) a 0 (44) 


ot ox oy ox od bx od ox oy oy 
= 3s os oT 
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ot u ox M oy wi oy i pd "by ne od ( ox i oy p Ox ¢ (2) 
Continuity 
an , 2 3 = 
: + ax (Ud) + Dy (Vd) = 0 (46) 
where 


U and V = depth-averaged horizontal velocity components at 
time t in the x- and y-directions, 
respectively, ft/sec 


n = displacement of the mean free surface with 
respect to the still-water level, ft 


p = mass density of seawater, silmaayizee 
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n + h = total water depth, ft 


€ and T = bottom friction stresses in the x- and 
bx by 2 
y-directions, respectively, lb/ft 
S 5 8 » and $ = radiation stresses which arise because of the 
xx xy Ws/ 


excess momentum flux due to waves (refer to 
Longuet-Higgins and Stewart (1964) for their 
significance), lb/ft 


T = lateral shear stress due to turbulent mixing, 
1b/£t2 


2S)/ 


a. CROSS-SECTION A-A 


OCEAN fe Breaker LINE 
BOUNDARY 


STILL- 
LINE / ee pee 
V. = 
A m A 
{ / 6 >o0h — a { 
[ D 
> (@) ap aD 
Sm 
4 
: : 
| \ 2 -*— SET-UP LINE 


b. PLAN 


Figure 8. Definition sketch for an irregular beach 
(swl = still-water level) 
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The condition yn > 0 is known as setup, and n <0 is called setdown. 
50. Bottom friction. At present, the numerical model uses a linear 


formulation for friction (Longuet-Higgins 1970). Thus 


Te 2oc <|(M |= U (47) 


Thy =Hoe <[uyepl? V (48) 


where c is a drag coefficient (of the order of 0.01) and <|u is the 


> 
eras 
time average, over one wave period, of the absolute value of the wave orbital 


velocity at the bottom. From linear wave theory 


2H 


<|u ~ WP gin Tan 


|> 


(49) 


orb 


Equations 47 and 48 are based on the assumption that the velocity components 
U and V of the current are small compared with the wave orbital velocity, 
<ju > 6 
| Be 
51. Radiation stresses. The radiation stresses are of major importance 
since they furnish the main forces for creating wave-induced currents. Refer- 


ring to Longuet-Higgins (1970), for monochromatic waves, they are defined in 


terms of the local wave climate as follows: 


Seeeene | (x -5) eos Q + (n - ) eine | (50) 


E n cos 9 sin 9 (51) 
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(n is the ratio of wave group celerity to phase celerity), 9 is the local 
wave direction (defined as shown in Figure 8), and E is the wave energy 
density. The values of H, k, and @ are obtained from RCPWAVE. 

52. Lateral shear. In the numerical model, the coordinate scheme is 
chosen such that x is positive in the offshore direction and y is approxi- 
mately in the alongshore direction. An eddy viscosity formulation is chosen 
for the lateral shear. The eddy viscosity is assumed to be anisotropic. 
Denoting - and ey as the eddy viscosities in x- and y-directions, respec— 
tively, in general, - is assumed to be a function of x and y and ¢ a 


y 
constant. Accordingly, 


= Bae & 
Txy o(c, ay fi x =) ee) 


For field applications, the eddy viscosity Ey is chosen according to the 


following relation given by Jonsson, Skovgaard, and Jacobsen (1974): 


2 
_, ele is 2 


€ 7 os 0 (55) 
4n h 


x 


This represents twice the value used by Thornton (1970). The value of ¢ 
was, in general, taken to be equal to the value of . at the deepest part 
(usually near the offshore boundary) of the numerical grid. 
Method of solution 

53. In view of the similarity among Equations 44-46 and the equations 
for long waves (Equations 7-9), CURRENT was developed by modifying WIFM. Thus 
CURRENT also is an implicit finite difference model and uses the SC scheme 
described previously. Details of the method of solution can be found in 
Vemulakonda (1984). 
Initial and boundary conditions 

54. In order to solve the problem under consideration, appropriate 
initial and boundary conditions must be specified. Usually an initial condi- 
tion of rest is chosen so that n » U , and V are zero at the start of the 
calculations. To avoid shock, the radiation stress gradients are gradually 
built up to their full values over a number of time-steps. The numerical 


computation is stopped when a steady state is deemed to have been reached. 
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55. The numerical model permits various types of boundary conditions 
among which are the following: 


a. "No flow" (wall). This type of boundary condition is used at 

FF closed boundaries such as the still-water line on beaches and 

at impermeable structures. The normal velocity is set to zero 
in this case. 


Io 


Uniform flux. In this type of open boundary condition, the 
flux at a boundary cell is made equal to that at the next 
interior cell. Thus the condition assumes 9(Ud)/ax = 0 or 
9(Vd)/d3y = 0 at the boundary. This type of condition is used 
for the lateral boundaries since it is a passive condition. 


c. Radiation. This open boundary condition requires that any 
transients developed initially inside the numerical grid should 
propagate out of the grid as gravity waves. It is of the form 
dn/ot + c(dn/dx) = 0 where c is the phase speed of a surface 
disturbance n(x,t) . It is often used by the wave-induced 
current model at the offshore boundary and is found preferable 
to a wall or constant elevation condition there. Both of the 
latter conditions are highly reflective, and, as a result, the 
transients tend to bounce back and forth between the offshore 
and nearshore boundaries and take a long time to damp out. On 
the other hand, the radiation condition seems to work quite 
well, allowing the transients to propagate out of the grid and 
permitting the setdown at the offshore boundary to assume an 
appropriate value. 

56. The boundary conditions frequently used in the wave-induced current 

model are illustrated in Figure 9. 

57. At present, the model allows for subgrid (thin-wal1) barriers such 
as jetties, provided they are impermeable and nonovertopping. The program 
essentially sets to zero the velocity component normal to the appropriate cell 


face. 
The Sediment Transport Model 


58. The sediment transport model predicts the transport, deposition, 
and erosion of noncohesive sediments such as sands in open coast areas as well 
as in the vicinity of tidal inlets. It accounts for both tides and wave ac- 
tion by using for input the results of WIFM, RCPWAVE, and CURRENT in terms of 
tidal elevations and currents, wave climate information, wave-induced cur- 
rents, and setups at the centers of grid cells. The model computes transport 
separately for straight open coast areas and areas in the vicinity of tidal 
inlets. In the case of the former, transports inside and outside the surf 


zone are treated separately. 
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Figure 9. Boundary conditions used in numerical model CURRENT 


Transport inside the surf zone 


59. Inside the surf zone it is the wave breaking process that is 
primarily responsible for the transport of sediment. This process is quite 
complex and not well understood. There is even considerable disagreement on 
the primary mode (bed load or suspended load) of sediment transport in the 
surf zone (Komar 1978). Thus a model that determines transport in the surf 
zone must be empirical, to some degree, in its formulation. 

60. The surf zone transport model used in this study is based upon an 
energetics concept developed by Bagnold (1963) who reasoned that the wave 
orbital motion provides a stress that moves sediment back and forth in an 
amount proportional to the local rate of energy dissipation. Although there 
is no net transport as a result of this motion, the sediment is in a dispersed 
and suspended state so that a steady current of arbitrary strength will trans— 
port the sediment. Thus breaking waves provide the power to support sand in a 
dispersed state (bed and suspended load), while a superimposed current (litto- 
ral, rip, tidal) produces net sand transport. 

61. The total littoral transport rate I, (vertically integrated and 
parallel to the shoreline) within the surf zone can be related to the wave 


conditions at the breaker line by 
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I = ie) sin cos a, (56) 
b 
where 
I = immersed weight sand transport rate (1b/sec) 
K = empirical coefficient 
On = breaker angle 


and the subscript b is used to denote conditions at the breaker line. 
62. Following Komar (1977), the local (vertically integrated) immersed 
weight longshore transport rate, per unit width in the cross-shore direction, 


may be written as 


tk, 2 
a Mee (0.5f£) ogy hv, (57) 
where 
ky = coefficient to be determined 
f = drag coefficient 
y= : = breaker index 
Moons local longshore velocity 
63. By integrating i, across the width of the surf zone Xp 
*b 
I) = f i, dx (58) 
0 
or 
“b 
tk 2 
< 1 H 
I, Sea, (0.5f£) oo J h v, dx (59) 


under the assumption that the coefficients ky and f are constants for a 


particular field site. Since the values of H , ) >» and h are known, 
being input to the sediment model, the integral on the right side of Equa- 
tion 59 may be determined numerically. For example, using the trapezoidal 


rule, 
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where IMAX corresponds to the number of water cells within the surf zone. 
Equation 60 allows for a gradual variation in cell size Ax . The velocity 
Vo is taken as the magnitude of the resultant of the total velocity 


components u, and v,, in the x- and y-directions. Thus 


T T 
Viens uw, + vs (61) 
where 
up =u +U (62) 
AE cals +V (63) 


For each computational grid line from the shoreline to the breaker line for 
each time-step, the value of I, is used to determine the unknown coefficient 


k, from Equations 56 and 59: 


2 
KH “e, sin Oh cos os 


iS ie we I, () 


The value of K is taken to be 0.39 if significant wave heights are used as 
in this study (Shore Protection Manual (SPM) 1984). 

64. Once kK is known, the local transport rate i, may be determined 
from Equation 57 and hence the local volumetric sediment transport rate dy > 


as in the following equation: 


2 
mk, foy hv, 


> ey oy & se) 


where 
0 = mass density of solids 


2.65 9 for sand 


ir 


jah) 
u 


ratio of volume of solids to total volume of sediment 


0.6 for sand 
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It is not necessary to know the value of f in order to solve for de in the 
above procedure. Once de is known, the local volumetric sediment transport 
rates q, and Ge for the cell may be determined by multiplying dp by 
Up/V» and ValVo >» respectively. 

Transport beyond the surf zone 

65. Beyond the surf zone, waves are not breaking. Currents (tidal, 
littoral, and rip) still transport sediment, but the sediment load is much 
smaller than the load in the surf zone. Waves still assist in providing power 
to support sand in a dispersed state. However, there is little turbulent en- 
ergy dissipation, and frictional energy dissipated on the bottom represents 
most of the energy dissipation. Bed load is the primary mode of sediment 
transport beyond the surf zone according to Thornton (1972). 

66. Since beyond the surf zone it is the tractive forces of currents 
(including wave orbital velocity currents) that produce sediment movement, a 
sediment transport by currents approach is taken. Again, since the complete 
physics of the problem is not completely understood, a semiempirical approach 
must be taken. In this model, the approach of Ackers and White (1973) is 
followed after appropriate modification for the influence of waves. 

67. Ackers and White (1973) studied sediment transport due to currents. 
They used the results of 925 individual sediment transport experiments to 
establish various empirical coefficients. The approach considers both 
suspended load and bed load. It is assumed that the rate of suspended load 
transport is dependent upon the total shear on the bed. Therefore, the shear 
velocity v, is the important velocity for suspended load transport. Bed 
load transport, however, is assumed to depend upon the actual shear stress on 
individual sediment grains. Ackers and White (1973) assume that this stress 
is comparable with the shear stress that would occur on a plane granular sur- 
face bed with the same mean stream velocity. Thus the mean velocity of flow 
v is the important velocity for bedload transport. 

68. Considering only currents (not waves), Ackers and White (1973) 
derived sediment transport rate in a dimensionless form. For convenience in 


practical application, this may be written as 


pele) Ses (66) 
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where 


q = total volumetric sediment transport rate per unit width normal to 
the current (vertically integrated combined bed and suspended 


sediment load CEES ECOTIES) 
Pp = porosity of sediment = 1 - a' 


D = sediment diameter which is exceeded in size by 65 percent (by 
weight) of the total sample 


n, = 1.0 - 0.2432 In Y (67) 


1/3 
Y=D [ees (68) 


v 
s = specific gravity of sediment 


y = kinematic viscosity of fluid 


C = exp [2-86 In Y - 0.4343 (In x)? - 8.128] (69) 
fra 22s n@ant 1) 
1/2: 
Y¢ 
_ 9.66 
m= + 1.34 (71) 
okval n-l 
v() (22 10s 224) 
Fe Ne (72) 


Re om? 


Equations 67, 69, 70, and 71 apply for 1 < Y < 60 (transition sediments). 


For values of Y greater than 60 (coarse sediments), C, n, , m, , and A 


1 1 


have the values of 0.025, 0, 1.5 and 0.17, respectively. 

69. Beyond the surf zone, both currents and nonbreaking waves exist. 
So the Ackers and White formulation derived originally for currents only must 
be modified for the presence of waves. The waves do not increase the level of 
turbulence since turbulence is confined to a narrow boundary layer by the 
oscillating wave orbital velocities. Since the shear velocity is dependent 
upon the intensity of turbulence and thus the total energy degradation rather 
than the net traction on individual sediment grains, the shear velocity is not 
changed by wave action. With the wave-induced turbulence confined to a narrow 


boundary layer and the waves propagating essentially without energy loss, the 
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effect of waves is to increase the traction on individual grains by increasing 
the mean velocity felt by the grains. Thus the mean velocity of flow must be 
increased, but the shear velocity must remain unchanged. The mean velocity of 
flow is increased by using the following equation developed by Bijker (1967) 
and modified by Swart (1974): 


a 2 
(v) = (v) i¢eig = (73) 
wave and current current 2 2; 
where 
2) 1/2 
Eo = Cy ea (74) 
x 10h 
Cc, = 18 log (4 ) (75) 
fw. = Jonsson's (1966) friction factor with D as bed roughness 
us = wave orbital velocity 


= <|| |> 


u 
orb 


Thus Equation 66 becomes 
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Equations 76 and 77 are used for calculating sediment transport beyond the 


surf zone. In these equations, v is interpreted as the total velocity Vo 


due to currents = te + v5 and Vi is obtained from the relation 
2 
Ps BVe 
%=V es [> (78) 
(e) C2 
VA 


where Pe is the bed shear stress and Cc. is the Chezy coefficient. From 
q , the local transport rates qe and a are obtained as before. 
Transport in the vicinity of inlets 

70. The flow and sediment transport in the vicinity of tidal inlets 
differ markedly from the flow and sediment transport in the surf zone for a 
straight open coast. The bathymetry in the inlet area is highly irregular 
with the presence of channels, bars, and shoals. The breaker line is gener- 
ally shifted farther offshore and is irregular. Breaking and decay of waves 
and wave-induced currents are the major mechanisms for transport of sediment 
in the surf zone near straight coasts, with tidal currents being of secondary 
importance. Generally Up is much less than Vr: In the vicinity of 
inlets, tidal currents are a major mechanism comparable to wave-induced cur- 
rents. Moreover, Un and Vp may be comparable. We are primarily inter- 
ested in the transport and deposition of sediment in the navigation channel. 
There is no guidance in the open literature as to how sediment transport in 
this area should be handled. In view of the factors mentioned previously, the 
model uses the Ackers and White formulation modified for the presence of waves 
(Equations 76 and 77) in this area. From previous experience (Vemulakonda 
et al. 1985), this approach was found to yield satisfactory results. 
Erosion and deposition 

71. In the case of noncohesive sediments, once the transport rates of 
sediment qd, and ay are known, changes in bed elevation can be determined 


from the continuity equation 


Si LS (79) 
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where & is the bed elevation. Equation 79 indicates that if more material 
enters a cell than leaves it, ¢ will increase (there will be deposition), 
and if more material leaves than enters, ¢ will decrease (there will be 
erosion). Equation 79 is applied in a finite difference form to all the grid 
cells at the end of each time-step to determine erosion and deposition. Note 
that an increase in ¢ means a decrease in still-water depth h and vice 
versa. Therefore, the values of h are updated simultaneously. 
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PART III: VERIFICATION AND BASE CONDITION TESTS 


Tides 


72. Astronomical tides are the primary driving force for currents 
within St. Marys Inlet; they also contribute significantly to the ocean cur- 
rents in the study area. WIFM is used to compute the currents for an average 
tide range in order to supply the sediment transport model with a time-series 
of depth-averaged horizontal velocity fields covering one tidal cycle 
(12.42 hr for the semidiurnal tide at Kings Bay). 

Verification 

73. Bathymetry. Most of the bathymetry and topography information used 
to define the grid cell elevations in Grid 1 (Figure 3) came from NOS nautical 
charts 11488, 11502, and 11503. Detailed soundings taken by the US Army Engi- 
neer District, Jacksonville (CESAJ), in June 1982 provided bathymetry for the 
navigation channels. All depths in the grid were referenced to mlw, and a 
datum difference of 3.0 ft between mlw and mean sea level (msl) was used. The 
maximum water depth in Grid 1 was 66 ft mlw. 

74. Prototype data. The prototype tide data used to calibrate and 
verify WIFM in this study consisted of tidal elevations and currents. Fig- 
ure 10 shows the locations of the tide and velocity gages deployed in the 
Kings Bay study area. Tide data were collected by the United States Geologi- 
cal Survey (USGS) and CEWES between September and December 1982. Currents 
were measured along ranges 1-4 (Figure 10) on 10 November 1982, and along 
ranges 5-7 on 12 November 1982. These surveys recorded approximately one 
tidal cycle. At each range, currents were measured at three stations: A, B, 
and C. At each station, velocities were measured at the surface, middepth, 
and close to the bottom. Only ranges 1-4 lay within the bounds of the compu- 
tational grids, so ranges 5-7 were not used in this report. These current 
measurements were accurate and error-free, so they were used by WIFM in veri- 
fication. The details of the prototype tide and current data collection 
effort are reported by Granat et al. (in preparation). 

75. Plates 1-4 show the measured tides at Gages 1-4 of Figure 10 for 
November 1982. The mean long-term tide ranges, as given by the 1982 NOS Tide 
Tables, vary between 5.8 ft (St. Marys Entrance, north jetty) and 6.0 ft 


(Fernandina Beach, Amelia River). The measured tide data for the range survey 
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Figure 10. Locations of field gages 
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date of 10 November 1982 agreed with these mean ranges and so represented an 
average tide for the study area. Since the prototype tides measured on 

10 November 1982 represented an average tide range, these tidal elevation sig- 
nals were used as boundary conditions in WIFM. The prototype range current 
data were used to verify the velocity computations. 

76. Plate 5 shows the prototype tide records for 10 November 1982. The 
sampling rate for the records was 5 min, and these data were spline filtered 
to remove high-frequency noise. Tides measured at Gages 2, 3, and 4 served as 
boundary conditions to the Amelia River, St. Marys River, and Cumberland Sound 
boundaries of the model. The signal from Gage 1, located at the south jetty 
of the inlet (Figure 10), was used as the boundary condition at the eastern 
edge of the computational grids. However, the travel time for a gravity wave 
between the eastern boundary and the actual location of Gage 1 is 25 min, so 
the boundary condition was phase shifted 25 min to account for this distance. 
The lateral ocean boundary conditions were interpolated between this offshore 
signal and the tide signal at the inlet (Gage 1). The boundary condition at 
St. Marys River (Gage 2) was also phase shifted 7 min to account for gravity 
wave travel time between the mouth of the river (the grid boundary) and the 
site of the prototype gage farther upstream. 

77. The zero datum shown in Plate 5 represents the mean for each 
measured tide record rather than a geophysical datum such as the National 
Geodetic Vertical Datum (NGVD) of 1929 . The elevations of the tide recorders 
used in this study were not referenced to a benchmark, so the relationships 
between the gage means are not known. The lack of a common datum caused 
numerous problems during calibration, since WIFM requires all elevations to be 
measured from a common datum. Since the tide gages were all fairly close to 
one another (less than 2 miles apart), even minor changes in elevations caused 
gradients great enough to change the flow patterns within the study area. 
These elevation adjustments were determined during the model calibration. 

78. Permeability of jetties. Since both jetties at St. Marys Entrance 
are awash at high tide and known to be permeable, the tidal model has to prop- 
erly simulate this effect on the velocity patterns. From field measurements 
taken by Florida Coastal Engineers in 1975, "it was estimated that up to 
28 percent of the total [tidal] flood flow enters [the inlet] through the per- 
meable jetties rather than at the ocean terminus of the structures." (Parchure 


1982, p. 27). Since the widths of the jetties are small compared with grid 
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cell dimensions, they can be modeled in WIFM as flow barriers placed at grid 
cell faces. The hydrodynamics of flow over these barriers is computed by the 
broad-crested weir formula (Chow 1959). The parameters of barrier submergence 
(head across the weir) and Manning's n in the formula dictate the flow rate 
or "permeability" across a barrier in WIFM. The permeable jetties at St. Marys 
Entrance were therefore simulated with submerged barriers in the tidal current 
model. 

79. An ad hoc method determined barrier "permeability" parameters for 
WIFM. Two initial assumptions were made to reduce the number of variables 
involved in parameter estimations. First, the crest of the submerged barriers 
used in WIFM was arbitrarily set to -4 ft msl. This depth ensured that the 
barriers would not become exposed during low tide. Second, it was assumed 
that the bottom friction in the study area, below -10 ft msl, could be 
approximated with a set Manning's n of 0.025. These assumptions reduce the 
variables affecting permeability to: (a) water velocity over the barrier, 
(b) water depth surrounding the barrier, and (c) the Manning's n of the bar- 
rier. The relationships between these variables were determined by a simple 
computational experiment. 

80. A horizontal flume with length scales of the same magnitude as 
St. Marys Inlet was modeled by WIFM. The flume is 16,000 ft long and 2,400 ft 
wide, and it has a submerged barrier obstructing half the channel width at the 
center of the flume. Plate 6 illustrates the plan view of the layout, and the 
velocity pattern for a typical computation. WIFM was run for 128 different 
combinations of flow velocity (1, 2, 3, and 4 fps), water depth (10, 20, and 
30 ft), and barrier Manning's n (varied from 0.025 to 0.050). Discharges 
per unit width were measured at the inflow and over the barrier for each run, 
and the permeability for the given conditions was computed as the percentage 
ratio of the latter to the former. It was determined that permeability was 
not a function of the flow velocity. 

81. Figure 11 shows the family of curves plotted from this experiment. 
To set a desired permeability for a jetty barrier in the tidal current model, 
the water depth at the jetty section is noted, and the appropriate Manning's 
n is determined by interpolating between isobath curves in Figure 11. The 
Manning's n values needed to simulate a 28 percent jetty permeability were 


determined for each barrier segment in this fashion. 
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Figure 11. The relationship of barrier permeability to depth and 
Manning's n 


82. Calibration of the tidal model required the adjustment of WIFM 
boundary conditions until the computed elevations at tide Gage 1 matched the 
prototype data for 10 November 1982. The model was then verified for 
correctness by successfully reproducing the velocities measured at ranges l, 
3, and 4. The WIFM boundary conditions were adjusted during calibration by 
adjusting the datums for the prototype tide signals and accounting for the 
phase differences in the signals due to their placements in the grid as input 
conditions. All of the datum adjustments were less than 2 in. Note that WIFM 
used a time-step of 60 sec for all the computations. 

83. Plate 7 shows where numerical gages were placed in Grids 1 and 2 in 
order to measure the computed velocities for the base and plane conditions in 
St. Marys Entrance. The gage sites in Plate 7 all correspond to either loca- 
tions where prototype data were collected during the survey of 10 November 
1982, or to important locations in the navigation channel. Table 1 equates 


the gage numbers in Plate 7 with the gage names used subsequently. The 
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missing numbers in the sequence correspond to gages outside the grid segment 
shown in Plate 7. 

84. Plate 8 shows the match between the computed and prototype tide for 
10 November 1982 at the four gage locations. The computed signals at 
Gages 2-4 are merely the prototype tide with the datum adjustments added, 
since these gages are boundary conditions in the model. Computations for tide 
Gage 1, in the inlet, match the prototype data. 

85. Plates 9-13 compare the model computations of tidal currents to the 
prototype surface and middepth velocities at ranges 1, 3, and 4. (Solid 
curves represent numerical results and dashed curves prototype data.) Varia- 
tions with time of both velocity magnitude and phase are shown over a tidal 
cycle. Since the numerical model is depth-averaged, in general its results 
would match the middepth measurements more closely. The agreement between the 
computations and the prototype data at the inlet (range 1) is excellent, both 
in magnitude and phase. The ability of the tidal current model to simulate 
the inlet velocities is crucial to the other aspects of this study, and the 
model performs this task well. In the case of ranges 3 and 4 (Plates 10-13) 
the numerical results represent the whole range. The computed and prototype 
velocities at range 4 (Cumberland Sound) also agree well. The velocity com- 
parisons at range 3 (St. Marys River) agree in magnitude but differ slightly 
in phase. This phenomenon is probably due to the drainage characteristics of 
large marsh areas around St. Marys River which lie outside the boundaries of 
the tidal model. 

86. Plates 14 and 15 show the computed velocity patterns in St. Marys 
Entrance for the peak flows of ebb and flood tide. The dashed portions of the 
barriers represent the permeable sections of the inlet jetties. The flow 
across the jetties can be seen on these Plates, and it appears that the flow 
is more pronounced across the south jetty. 

87. In summary, the tide model used prototype gage elevation data for 
forcing boundary conditions. The measured tidal elevation at the south jetty 
of the inlet was reproduced in the numerical model. There was good agreement 
of numerical results with measured velocity data at range 1 in the inlet and 
satisfactory agreement at interior velocity ranges. Where the flows are 
influenced by other features in the region interior to the inlet, such as 
marshes which are not included in the tide model of Model B, close agreement 


is not expected. This lack of agreement should not cause concern since the 
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interior flows are studied by Model A and since the main purpose of Model B is 
to study coastal processes in the region mostly exterior to the inlet. There- 
fore, the calibration and verification of the tide model are complete and 
successful. 

88. Since the channel bathymetry and geometry used in verification 
tests are close to existing (base) condition and the tide of 10 November 1982 
is representative of the mean tide, the results of WIFM from verification 
tests were also taken to be those for base condition. They were used accord- 
ingly in the sediment transport model. The reader should note that the tidal 
conditions of 10 November 1982 were used in Model A also for base condition. 
The tide model generated a data file, consisting of tidal elevations and 
velocities, for each grid cell for each half hour of an approximated semi- 


diurnal period of 12.50 hr for later use in the sediment model. 


Waves and Wave-Induced Currents 


89. The hydrodynamic models RCPWAVE and CURRENT were extensively tested 
and compared with analytic solutions, laboratory data, and available field 
data during their development (Ebersole, Cialone, and Prater 1986 and 
Vemulakonda 1984). Considerable experience has been gained previously at WES 
in field application of these models (Vemulakonda et al. 1985). So reliance 
can be placed on the results of these models. The models do not require site- 
specific calibration. Because synoptic field data on waves and wave—induced 
currents were unavailable for the project area, no separate verification tests 
were performed for these models except indirectly through sediment model veri- 
fication. The models were run for the base condition using the same bathyme- 
try and channel geometry as in the tidal model. The results of the models 
were used in verification and base tests of the sediment transport model. 

Wave climate 

90. One of the primary objectives of the wave and wave-induced current 
model runs is to furnish input to the sediment transport model. In the case 
of the sediment transport model, the interest is in sediment transport and 
yearly shoaling rates in the navigation channel under an average year's wave 
climate, including normal storms but excluding extreme storms such as hurri- 
canes and other tropical storms. So the wave climate for an average year 


at the project site was obtained from the WES Wave Information Study (WESWIS) 
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based on 20-year hindcast. This information was in the form of frequency of 
occurrence of waves in terms of predominant direction, significant wave 
height, and period bands in a depth of 60 ft mlw. Table 2 shows a sample of 
WESWIS data for St. Marys Inlet. The wave approach angle notation in this 
table is different from that used in the rest of this report. Angles in the 
table are measured with respect to the shoreline. Consider waves with an ap- 
proach angle of 70.0 to 79.9 deg and significant wave heights in the band 0.0 
to 0.49 m. They are distributed in period bands between 0.0 to 11.0 sec and 
greater. The total frequency of occurrence of these waves summed over all the 
period bands is 3.515 percent (3,515 + 1,000) or 0.03515. Similarly, WESWIS 
provides wave information in direction bands of 10 deg from 0 to 180 deg for 
all the wave height bands (0 to 5.00 m and greater). 

91. In this study, these data were further consolidated into 79 differ- 
ent incident wave conditions (combinations of significant wave height, period, 
and direction) to run the wave and wave-induced current models. For conve- 
nience in running the sediment transport model, wave condition 80 was defined 
as a null wave condition when there was no significant wave activity. These 
combinations are listed in Table 3 which shows the percentage of occurrence of 
each condition. The directions represent angles in degrees measured from azi- 
muth 87.5 deg (approximate shore normal direction). Negative angles signify 
waves coming from directions south of the normal; positive anglés signify 
waves coming from directions north of the normal. The wave combinations shown 
in Table 3 are obtained from Table 2. Consider the example from Table 2 
again. Since the wave approach angle is between 70.0 and 79.9 deg, the aver- 
age value of 75.0 deg is taken. In terms of the notation of Table 3, the wave 
direction becomes -12.5 deg. Since the wave height band is from 0.0 to 
0.49 m, the mean value of 0.25 m or 0.82 ft is taken for the significant wave 
height. As for period, on the basis of the distribution of Table 2, a mean 
period of 8.0 sec is taken. These are the values shown for wave 28 in 
Table 3. 

Jetties 

92. To account for diffraction of waves due to the two jetties of 
St. Marys Inlet, a special subroutine was developed. It used the diffraction 
solution of Penney and Price (1952). The wave model was first run without 
accounting for the presence of the jetties. The diffraction subroutine took 


the solution near the jetties as input and modified it to allow for 
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diffraction around the jetties. For this, the actual layout of the jetties is 
used. The procedure was somewhat similar to that of Perlin and Dean (1983). 
During the development of the subroutine and the procedure, several tests were 
performed including comparison of its results to the laboratory data of Hales 
(1980) for a single structure case, and to two physical hydraulic model tests 
conducted at CEWES for the two jetty case. In each case, the results of the 
subroutine compared favorably with laboratory data. In the grid for CURRENT, 
the jetties were represented in a stair-step fashion similar to that in WIFM. 
CURRENT treated them as thin-walled nonovertopping impermeable barriers. 

93. Because of the highly variable nature of the computational grid, 
the wave model was run on a uniform grid with 500-ft by 500-ft cells, and its 
results were interpolated to the variable grid. The wave and wave-induced 
current models were run for each of the 79 wave conditions. There are no 
waves or wave-induced currents corresponding to wave 80. Each of the wave 
conditions represented the offshore boundary condition for the wave model. 

The model was run for the condition, and its results were stored in the form 
of wave height, direction, and wave number at each grid cell. They were next 
used as input to CURRENT which computed and stored on a file the setup nN and 
the two velocity components U and V for each grid cell for each wave. For 
convenience the corresponding wave information for each cell was also stored 
on the same file. Note that in general CURRENT used a time-step of 50 sec and 
in each run calculations were continued until an approximate steady-state con- 
dition was reached by the current field. 

Results 

94. For convenience, results for only three typical cases out of the 79 
listed in Table 3 will be presented here. They have been selected so that 
they represent waves coming from south and north of the shore-normal direction 
and approximately along the shore-normal direction. It is convenient to pre- 
sent the results from the wave and wave-induced current models in terms of the 
uniform grid in the computational plane rather than the variable grid (Fig- 
ure 3). One advantage of this type of display is that the results for the 
entire grid can be shown on an 8-1/2-in. by 1l-in. sheet of paper. However, 
there is a disadvantage in that cell centers are not at the proper distances 
relative to each other. Thus, boundary cells appear much closer to the center 
than they really are. Moreover, the cell dimensions are distorted. Cells 


close to the inlet, the barrier islands, and the navigation channel appear to 
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be relatively larger; and as one moves away from this region (for example, 
closer to the lateral and offshore boundaries) the cells appear to be rela- 
tively smaller than they really are. In what follows, for convenience, the 
results will be shown on the uniform grid. 
95. Figure 12 shows the region covered by the 50- by 73-cell uniform 
grid in the computational plane. The grid is 50 cells wide in the alongshore 
75 
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Figure 12. Uniform grid and bathymetry in 
computational plane 

direction and 73 cells long in the onshore/offshore direction. The figure 
shows Amelia and Cumberland Islands, the navigation channel, the locations of 
the two jetties on St. Marys Inlet (the jetties are stair-stepped for the CUR- 
RENT model), and bathymetric contours with elevations referenced to msl (mlw 
plus 3 ft). The offshore boundary of the grid is at an approximate depth of 
63 ft msl. Note the shoals offshore, south of the south jetty, and north of 
the north jetty. For convenience, these will be referred to hereafter as the 


offshore, south, and north shoals, respectively. 
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96. Figures 13, 15, and 17 display the results of the wave model, at 
cell centers, corresponding to waves coming from three different directions. 
These are waves 22, 45, and 59 from Table 3. These three conditions will be 
referred to as cases A, B, and C, respectively. For all three cases, the 
significant wave height in 63-ft depth of water msl is identical and equal to' 
7.4 £t. The period is roughly the same. The wave directions are quite 
different. In each of the figures, the length of an arrow (vector) is 
proportional to the wave height (a scale is shown), and the direction of the 
arrow indicates the direction in which the waves are progressing. For 
clarity, only vectors for alternate cells in each coordinate direction are 
plotted. 

97. Figures 14, 16, and 18 present the wave-induced currents at grid 
cell centers corresponding to cases A, B, and C. In these figures, the length 
of an arrow is proportional to the magnitude of the current (a scale is shown), 
and the direction of the arrow indicates the direction of the current. The 
currents are depth-averaged. For convenience and clarity, only vectors in 
alternate cells in each coordinate direction are plotted. 

98. Figure 13 corresponds to a wave of period 7.2 sec coming from 
azimuth 120 deg in 63-ft depth msl (Case A). The waves respond to the off- 
shore shoal. The wave height increases, and the wave direction changes as the 
waves go over the shoal. The wave height decreases, and the waves resume 
their original direction once the waves pass the shoal. The waves converge on 
the south shoal due to refraction, move parallel to the jetty, and break on 
the shoal. Because of the sheltering effect of the south jetty, very little 
of the incident wave energy goes past the jetty tips into the inlet. Note 
also the sheltering effect behind the north jetty resulting in very little 
wave action there. The waves converge on the north shoal, and the wave energy 
spreads out (diverges) due to a "bay" effect as the waves reach the shoreline 
of Cumberland Island. Near the approximately straight shorelines of both bar- 
rier islands, the wave height decreases because of wave breaking and decay. 

99. Figure 14 shows the wave-induced currents corresponding to Case A. 
Near the straight part of the shorelines of Amelia and Cumberland Islands, the 
currents are mainly parallel to the shore and move to the north. However, 
near the south shoal, because of wave refraction and breaking, the currents 
tend to move in a westerly direction. The net result is the counterclockwise 


circulation we see over the shoal. The currents are the largest in this 
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Figure 14. Wave-induced 
currents for Case A 
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Figure 16. Wave-induced 
currents for Case B 
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region. Currents are strong on the inside of the north jetty because waves 
advance and break along the jetty. These currents have a westerly direction 
and advance into the inlet. Because of diffraction, currents are very weak 
behind the north jetty. 

100. Figure 15 corresponds to a wave of period 7.8 sec coming from 
azimuth 80 deg in 63-ft depth of water msl (Case B). In this case, since the 
waves are approximately normal to the shoreline and the offshore contours, 
there is not much refraction of the waves offshore or even near the straight 
line portion of the shoreline. The waves converge on the south shoal, because 
of refraction, resulting in higher wave heights on the shoal. There is a 
similar convergence on the north shoal and a small divergence of wave energy 
near Cumberland Island. The incident wave direction is such that there is 
very little sheltering due to the two jetties. As a result the waves propa- 
gate straight and far into the inlet because the depth contours are approxi- 
mately straight and parallel to the waves inside the jetties. The wave 
heights are large between the jetties. 

101. Figure 16 displays the wave-induced currents for Case B. In this 
case, because the incident waves are approximately normal to the shoreline, 
there are no noticeable currents along the straight portions of the shoreline. 
Because of wave convergence and breaking, the currents are strong over the 
south shoal. A circulation pattern may be observed on the shoal. As the 
waves propagate straight and unchanged between the jetties without breaking or 
decaying, there are no noticeable wave-induced currents in this region. Cur- 
rents may be observed on the north shoal because of wave convergence, break- 
ing, and decay there. These currents are smaller than those observed on the 
south shoal. 

102. Figure 17 corresponds to a wave of period 6.9 sec and azimuth 
60 deg in 63-ft depth of water msl (Case C). The waves refract on the off- 
shore shoal. They refract and converge strongly on the north and south 
shoals, resulting in higher wave heights on both shoals. Since the waves are 
aligned approximately parallel to the two jetties, there is very little 
sheltering due to the jetties so that the waves propagate deep into the area 
between the jetties. They break and decay near the straight portions of the 
shoreline. 

103. Figure 18 represents the wave-induced currents for Case C. Near 


the straight reaches of the shoreline the currents are parallel to the shore 
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and in the southerly direction, as one would expect. The currents are strong 
over the north and south shoals because of wave breaking and decay. The pat- 
tern of the currents is complicated. Currents move in an easterly direction 
along the interior of the north jetty and westerly direction along the 
interior of the south jetty. 

104. In summary, the overall results of the wave and wave-induced cur- 
rent models used for verification and base conditions are reasonable and 
behave in a manner one would expect, given the complicated bathymetry of 
St. Marys Inlet region and the two jetties on the inlet. The incident waves 
respond differently to the bathymetry, the shoals and the jetties, depending 
on their direction of incidence. The wave-induced currents depend on the 
bathymetry, the waves everywhere in the grid, and whether or not the waves 


break and decay in a given region of the grid. 


Sediment Transport 


Verification tests 

105. In order to make a strict verification of the sediment transport 
model, it is necessary to have either long-term (several years long) informa- 
tion on shoaling rates in the navigation channel and bathymetric changes in 
the general area or actual wave measurements made simultaneously with measure- 
ments on shoaling rates and bathymetric changes over a shorter time period (a 
few months). The latter type of data are not available for the project area. 
As for the former, examination surveys are available for the channel. As men- 
tioned previously, the approach used by the sediment transport model does not 
account for extreme storms. So the prototype data selected should not include 
periods of such storms. As for dredging, it is possible to simulate dredging 
in the numerical model provided detailed information is available on the loca- 
tions and durations of dredging and the amounts of material dredged at each 
location. Usually, such detailed information in terms of computational grid 
cells is not available from dredging records. Therefore the prototype data 
should not include periods of dredging. In the case of St. Marys Inlet, the 
navigation channel was deepened to the existing condition (40-ft project 
depth) in 1978 and 1979 so only 5 to 6 years of prototype data are available. 
The channel still has not stabilized after the deepening. Sediment transport 


and other processes continue to be in a state of transition. Out of the 
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available information on examination surveys for the navigation channel, we 
were able to locate only one set of examination surveys covering a period of 
approximately 1 year (1980 through 1981) which was free from the effects of 
dredging and severe storms. The duration of this data set is too short for 
the data set to be used for strict verification of Model B results which are 
based on 20-year hindcast wave data. Therefore, a strict verification of 
Model B results with the data set is not possible. Instead, the average 
yearly erosion/deposition rates along the channel obtained from the data set 
will be compared with Model B results to see if the numerical model results 
are reasonable and agree with the trends and shoaling magnitudes exhibited by 
the field data. 

106. Prototype data. The field data set consisted of seven examination 
surveys conducted by CESAJ during 1980 and 1981 between sta -80+00 and 
sta 325+00 (Figure 19). For convenience, this pre-1985 CESAJ stationing will 


be used throughout this report. The locations and dates of the surveys are 
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Figure 19. Pre-1985 CESAJ stationing 
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shown in Table 4. On the basis of several tests, it was determined the datum 
used in survey 2 was in error by 0.5 ft. This is not surprising since the 
reach of channel surveyed was far away from the tide gages used to locate the 
datum. The datum for this survey was adjusted accordingly. 

107. The field data were examined in two ways. First, surveys l, 2, 
and 7 were used to determine average yearly erosion/deposition rates. At each 
of the 35 locations along the channel corresponding to computational grid cell 
centers, depths across the width of the channel were averaged, and the 
erosion/deposition rates were computed and extrapolated to feet/year values. 
(Following a similar procedure, but computing the average depth for each cell 
from 16 spatially distributed points in the cell, yielded results that were 
close to the results obtained from averaging the cross-section depths). Next, 
the total period was broken down into three separate periods of approximately 
4, 2, and 6 months, based on the survey dates. At each of the 35 locations, 
the erosion/deposition rates obtained for these periods were converted to 
feet/year values, and the extreme values at each location were determined. 


Figure 20 is a plot of the average and extreme values from the prototype data 
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Figure 20. Prototype data on erosion/deposition rates 
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at different stations along the channel. The sign convention that erosion 
rates are positive and deposition rates are negative is used hereafter. 

108. Testing procedure. At the start, the sediment model used the 
bathymetric information from field surveys 1 and 2 for the channel. Outside 
the channel, the bathymetry used was identical to that used by the tide, wave, 
and wave-induced current models for base condition since better detailed 
information was not available. 

109. To generate a wave sequence for 1 year for the verification and 
base tests of the sediment model from the waves given in Table 3, the follow- 
ing procedure was used. Each wave event in the sequence was assumed to be 
steady with a duration of 4 hr. This is a reasonable assumption from field 
experience and measurements, provided extreme storms are ruled out as done 
here. Each wave condition (1 to 80) of Table 3 was identified with a fre- 
quency of occurrence. During the running of the sediment model, wave condi- 
tions were selected such that each of the 80 conditions occurred at the 
frequency shown in Table 2. Thus, the waves used by the sediment model re- 
flected nature in terms of wave statistics provided by WESWIS. The same waves 
were used for base and Plan 1 tests. 

110. The sediment model used a time-step of 1 hr. This value was 
considered optimum on the basis of testing and previous experience. The 
computational sequence employed by the sediment model consisted of the fol- 
lowing steps: 

a. Read in the local bathymetry. 
b. Pick the first wave condition. 


c. Read in the corresponding wave information (wave height, 
angle, period, wave-induced velocities, and setups/setdowns). 


. Read in the first hour of the tide data (tidal velocities and 
elevations). 


- Combine the above quantities to obtain a total velocity field, 
wave field, and local depth. 


f£. Compute sediment transport quantities and the associated ero- 
sion and deposition rates. 


g- Repeat steps d, e, and f at 1l-hr intervals for a total of 
4 hr. 


h. Pick the next wave condition and continue steps c through g, 
Gie@o 


As indicated previously, the local still-water depth h for each cell was 


updated at each time-step based on the erosion or deposition in the cell. The 
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total local depth, which is the sum of h , ", and 1 , was also updated. 
The total velocity components Up and Vr were adjusted on the basis of sim- 
ple continuity to account for the change in bed elevation of the cell. It was 
observed from running the sediment transport model that model results in terms 
of erosion/deposition rates (ft/year) along the channel became approximately 
constant after the model was run for 150 to 180 prototype days. There were 
minor variations from run to run as the total time was increased, but the 
trends and magnitudes stabilized. Therefore, the above sequence was performed 
for 180 instead of 365 prototype days to compare with field data for 
verification. 

111. Results. Based on the trends exhibited by the prototype data, the 
reach of the channel between sta -—80+00 and sta 325+00 was divided into 
seven zones for verification (Figure 21). A similar approach was used for 
Model A verification. Four to six computational grid cells were in each zone. 
The value assigned to a zone is the average of the values for the cells in the 
zone. A comparison of the prototype average erosion/deposition rates with 
results of Model B is shown by zones in Figure 22. Also shown in the figure 
are prototype extrema based on 1 year of prototype data. Model B results show 
the same trends as the prototype average results and are in approximate 
quantitative agreement in zones 1-4 (between sta -80+00 and sta 241+56). It 
is not surprising that they do not match quite as well between sta 241+56 and 
sta 325+00. This is a highly dynamic region, especially outside the jetty 
tips (sta 251+00) and is very much dependent on the actual (rather than aver- 
age) wave climate that existed between surveys. This can be seen in the large 
spread between prototype extrema. There is movement of material from the off- 
shore bar and shoals into the channel. As was pointed out previously, proto- 
type data of 1 year's duration are not necessarily representative of a 20-year 
average. On the whole, Model B results are reasonable and in agreement with 
field data. 
Base tests 

112. The only difference in the bathymetries used at the start of veri- 
fication and base tests was in the navigation channel. The sediment model 
used the channel bathymetry from the CESAJ survey of June 1982 for the base 
tests because the same bathymetry was used in all the other models for base 


condition. 
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Figure 21. Zone numbers assigned to channel reaches for verification tests 
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Figure 22. Comparison of prototype data with Model B results 
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113. At the time Model B computations were made, up-to-date field sur- 
vey information was not available on channel bathymetry between sta 325+00 and 
sta 399+74, nor was up-to-date bathymetric information available for areas on 
either side of the channel. Model B used the best available information, 
which usually was CESAJ construction dredging survey information and the 
bathymetric information from NOS charts. Unfortunately, this information does 
not seem to represent the current bathymetry for the reach of the channel be- 
tween sta 325+00 and sta 399+74 and the areas of either side of the channel in 
this reach, according to the latest CESAJ surveys. These surveys seem to 
indicate that the depths may be greater (by 5 ft or more) in this reach. For 
this reach of the channel, the information available to us at the time of com- 
putations indicated the depths were greater than the existing project depth of 
40 ft and that the channel seemed to be in a state of erosion; therefore it 
did not require maintenance dredging. There was no quantitative information 
available on erosion/deposition rates for this reach. 

114. The sediment transport model followed the same testing procedure 
as it did for verification. The model was run for 200 prototype days, and its 
results were converted to channel erosion/deposition rates (ft/year). For 
base tests, the entire length of channel offshore of sta 399+74 was consid- 
ered. The channel was divided into ten zones (zones 1 to 10 in Figure 23). 
Figure 24 shows the erosion/deposition rates by zone for base. When the 
results for verification and base are compared (Figures 22 and 25), it is 
observed that there is similarity in the trends and magnitudes. This is not 


surprising since the channel bathymetry in the two cases is not that much 
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Figure 23. Zone numbers assigned to channel reaches for 
base and Plan 1 results 
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different inside the jetties and the forcing functions (tides, waves, and 
wave-induced currents) are the same for both cases. In both cases, there is 
deposition outside the jetty tips (sta 251+00). It changes to erosion inte- 
rior to the jetty tips because of circulation due to wave-induced currents. 
The heaviest deposition rates are observed near the jetty tips. This is the 
area where the channel cuts through the offshore bar and material from the bar 


tends to move into the channel and deposit there. 
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Figure 24. Erosion/deposition rates (ft/year) for base condition 
from Model B (+ = erosion, - = deposition) 
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PART IV: PLAN CONDITION TESTS 


Paltaniwel! 


115. Model B tested only one plan condition which will be referred to 
hereafter as "Plan 1." The plan is to (a) widen the navigation channel to 
500 ft, with the widening taking place on the north side of the present 
entrance and offshore channels; (b) extend the channel on the ocean side, with 
the extension being at an angle 20 deg south of the present channel center 
line at sta -97+76 approximately; and (c) deepen the channel to -49 ft mlw 
(46-ft project depth plus 3-ft advance maintenance). The channel is to have a 
trapezoidal cross section with side slopes of 3H:1V. Figure 26 shows details 
of the planned channel layout and cross section. As requested by the Officer 
In Charge of Construction (OICC), TRIDENT, the plan tested assumes also that 
the landward 1,000 ft of the south jetty will be made sand-tight 
simultaneously. 

116. In view of the urgent need expressed by OICC for Plan 1 results 
from Model B for design of the entrance and offshore channels, the wave and 
wave-induced current models were not rerun for the Plan 1 condition as origi- 
nally planned. Running the models again would have delayed the results con- 
siderably. Moreover, since the changes from base to Plan 1 condition of the 
navigation channel were reflected mainly in the cell size and bathymetry for 
one row of cells in the computational grid, it was felt the effect of channel 
modification on waves and wave-induced currents would be minor compared to its 


effect on tides and sediment transport. 


Computational Grid 


117. As indicated previously in Part II, the computational grid for 
Plan 1 (Grid 2) retained the major features (overall dimensions, orientation, 
number of cells, etc.) of the grid for base (Grid 1). The only difference 
between the two grids lies in the mapping of the row of cells corresponding to 
the navigation channel. These cells were made 500 ft wide by minor adjust- 
ments of the cells on either side. In view of the rectilinear nature of the 
grid, the navigation channel was represented in a stair-step fashion where it 


turned south. It was assumed that dredging for the navigation channel 
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extension stopped wherever the natural ocean depth became equal to or greater 
than 49 ft mlw. (This location of the oceanward entrance of the navigation 
channel would be determined in the field from the latest bathymetric surveys 
for the final channel design.) In the navigation channel itself, the planned 
channel depths were used. The bathymetry used outside the channel was the 


same as that for base condition. 


Tides 


118. To properly model the sand-tightened section of the south jetty in 
the numerical model, the crests of the barriers simulating this jetty section 
were raised to the prototype jetty elevation of +3 ft msl; and the Manning's 
n values governing flows over the barriers were changed appropriately. 

119. Plates 16-19 compare the computed base and Plan 1 velocities 
(magnitudes and phases) at seven sites in the inlet (refer to Table 1 and 
Plate 7 for locations of these sites). All of the changes in tidal currents 
are due to sand-tightening of the south jetty. The peak velocity at tide 
Gage 1 (Plate 16) has increased by approximately 10 percent between base and. 
Plan 1 due to sealing of a section of the south jetty. The gages at the 
throat of the inlet (Endeco velocity Gage 2, range survey Gages 1-A, 1-B, and 
1-C, and Fort Clinch) (Plates 16-19) show negligible change in velocity. The 
velocity at the ocean end of jetties (Plate 18) increases by about 10 percent 
in both ebb and flood for Plan 1 and shows a slight phase shift. 

120. Plates 20-21 show the tidal current patterns near the inlet for 
maximum ebb and flood for Plan 1. For clarity, the plotting of velocities 
below 0.1 fps is suppressed in these figures. These two plates can be com- 
pared to the base condition patterns of Plates 14-15, but few differences are 
apparent in a visual examination. For convenience the vector differences 
between Plan 1 and base condition velocities are shown in Plates 22-23 for the 
same region near the inlet, for maximum ebb and flood, respectively. Note the 
change in velocity scale. The plotting of velocity differences below 0.05 fps 
is suppressed in these figures. Both figures indicate that sealing the south 
jetty exerts local changes on tidal currents. The large difference vectors at 
the landward end of the south jetty represent a decrease in velocities between 
base and Plan 1 since flows in this area are stopped by sealing the jetty. No 


other significant changes in the current patterns were noted within the study 
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area for Plan 1. The extension of the navigation channel at the seaward end 


produced almost no effect upon the current patterns. 


Sediment Transport 


Testing procedure 
121. The testing procedure used was similar to that for base conditions 


except the computations were performed with Grid 2, and the bathymetry at the 
start of computations corresponded to Plan 1 conditions. The sediment trans- 
port model used the results of the tide model for Plan 1 and the results of 
the wave and wave-induced current models for base conditions. As for the base 
test, computations were performed for 200 days of prototype time, and the 
results were used to estimate yearly erosion/deposition rates along the 
channel. 
Results 

122. The channel was divided into 11 zones for Plan 1 (Figure 23). The 
exact offshore limit of zone 1A was yet to be determined from field surveys. 
Figure 27 shows the erosion/deposition rates by zone for Plan 1. A comparison 
of base and Plan 1 results is plotted in Figure 28. The model predicts an 
increase in both deposition rates and erosion rates between sta -—9/7+76 and 
sta 325+00 from base to Plan 1. Zone 1A is not shown in the figure. This 
zone indicates on the average a slight erosional tendency with rates of the 
order of 0.1 ft/ year or less. The model predicts deposition in zones 8 and 9 
(sta 323402 to sta 374+94) for both base and Plan 1. For Plan 1, the deposi- 
tion rates are of the order of 1.0 to 1.4 ft/year. Model B predicts large 
erosion rates in zone 10 (sta 374+94 to sta 399+74). Since there is no quan- 
titative field information on sedimentation rates in zones 8-10, it is diffi- 
cult to comment on Model B predictions for this reach. It is suspected that 
since the bathymetric information used for this channel reach and adjacent 
areas was not up-to-date, it might have caused deviation of Model B results 
from field experience. Another contributory factor might be the grain sizes 
found in this reach, which are much larger than elsewhere in the study area. 
The model assumes the same grain size distribution throughout the study area. 
The local deviation of grain sizes might have resulted in the prediction of 
larger erosion and deposition rates locally. There is reason to believe the 


effect of these factors is restricted to Model B predictions for zones 8-10 
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Figure 27. Erosion/deposition rates (ft/year) for Plan 1 condition 
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Figure 28. Comparison by zones of erosion/deposition rates 


for base and Plan 1 
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and does not extend to the model results for the rest of the study area. 

123. In terms of yearly shoaling volumes, if the reach of channel from 
sta -80+00 to sta 325+00 is considered, the results translate to approximately 
475,000 cu yd/year for base and 788,000 cu yd/year for Plan 1 allowing for the 
wider Plan 1 channel, or an increase of approximately 66 percent from base to 
Plan 1. The base volume is of the same order as the maintenance dredging vol- 


umes recorded in CESAJ dredging logs. 
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PART V: MODELING LIMITATIONS 


124. The numerical models used in this study were the most advanced 
models available at the time the study was undertaken. However, they do have 
certain limitations which must be kept in mind in order to view the results 
obtained from this study in proper perspective. As previously indicated, 
numerical models represent an approximation to the physical processes. The 
degree of approximation depends on the physics in the formulation of the indi- 
vidual models, the resolution of the computational grid, and the time-step 
used in computation. The assumptions made in the models and the limitations 
of the models have been given previously along with the description of the 
models in Part II. In this study, the computational grid resolution and the 
computational time-steps used have been chosen, on the basis of experience and 
testing, such that the results obtained would be reasonably accurate for engi- 
neering purposes and, yet, the computational costs would not be prohibitive. 

125. Generally, the hydrodynamic models are more exact than the sedi- 
ment transport model because more insight into the hydrodynamics is available 
and more experience has been gained in modeling the hydrodynamics numerically. 
With proper calibration and verification, tidal hydrodynamic models such as 
WIFM can predict tidal elevations very accurately and tidal currents fairly 
accurately. Monochromatic wave models such as RCPWAVE are fairly accurate in 
open coast areas. Their results near structures and inlets are more approxi- 
mate because of the difficulties and expense in modeling diffraction near 
structures and wave/current interaction near inlets. As for wave-induced cur- 
rent models such as CURRENT, people have less experience with them than with 
tidal and wave models. Wave-induced current models are reasonably validated 
for open-coast situations. Their results are more approximate near inlets and 
jetties because the hydrodynamic processes are more complicated and less 
understood, the wave fields are less accurately known, and there is a lack of 
field data to validate the models. 

126. Sediment transport is the most important aspect of the project for 
project design; yet the sediment transport model is the least exact of all the 
models, and the uncertainty is the greatest with this model. The uncertainty 
exists because sediment transport in general involves complex interactions 
between the bed and the flow which are not well understood. Of all types of 


sediment transport, sediment transport near inlets under the combined action 
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of waves and currents, as typified by this study, is one of the most compli- 
cated and least understood processes. The sediment transport model employed 
in this study uses fairly simple empirical formulas which are based on labora- 
tory and field data. It reflects the inaccuracies inherent in the formulas as 
well as the inaccuracies in the results of the three numerical hydrodynamic 
models. 

127. In this study, a mean tide, and an average year's wave climate 
based on 20-year averaging of wave statistics were used in running the sedi- 
ment transport model to estimate the yearly shoaling rates in the navigation 
channel. In reality, the tidal cycle is more complex, involving spring and 
neap tides; and the wave climate varies from day to day, season to season, and 
year to year. Severe storms such as hurricanes, which have a dramatic impact 
on sediment transport and channel shoaling, have been excluded from this 
study. As a result, sediment transport and channel shoaling rates in any 
given year may deviate significantly from the values predicted in this study. 
Moreover, short-term rates such as averages over a month or a season may dif- 
fer markedly from average rates over a year. Even the nature of sedimentation 
at a particular location may change from erosion to deposition and vice versa. 
This change is exemplified by the field data on shoaling rates shown in Fig- 
ure 20. Therefore, the results of this study will provide reasonable esti- 
mates of the long-term yearly average values of sediment transport and channel 
shoaling rates, provided severe storms are excluded and the uncertainty in the 
results for the reach of channel between sta 325+00 and sta 399+74 is noted. 

128. In general, the uncertainty in the predictions of the sediment 
transport model is reduced by verifying the model with field data from the 
project site. For St. Marys Inlet, field data from navigation channel surveys 
were available for about a year (November 1980 to December 1981) and are free 
from the effects of extreme storms. The model verification, as shown in Fig- 
ure 22, is good. In a sense, the verification shown is an indirect verifica- 
tion of the modeling system approach as a whole. The yearly shoaling volumes 
predicted by the model for existing conditions are comparable to yearly main- 
tenance dredging volumes recorded by CESAJ. In general, with proper verifica- 
tion numerical sediment transport models are better at predicting the effect 
of a change from one condition to another, such as from base to plan, than at 
predicting an absolute condition such as base or plan alone. In view of these 


facts, it is estimated that Model B results on sediment transport are accruate 


US 


to within +25 percent for base and Plan 1 conditions. 

129. To keep things in perspective, it should be pointed out that at 
present the only possible alternative to a numerical sediment transport model 
‘ls a physical movable-bed hydraulic model. Movable-bed coastal models are 
fairly complicated and expensive to construct and operate. Such models 
require more time than a numerical modeling effort. At present, there is no 
universal agreement on the scaling relations to be used. Movable-bed coastal 
models involving the combined action of waves and currents near inlets (the 
type required by the present study) are the most complicated of all coastal 
models and the least understood. Their results are approximate because there 
has to be a compromise between scaling relations necessary for waves only and 
scaling relations required for currents only. According to established ex- 
perts, the accuracy of a coastal movable bed model of this type in the hands 
of an expert will be of the order of +50-100 percent. Thus, the results of 
the numerical modeling system employed in this study are definitely better 


than the alternative. 
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PART VI: RECOMMENDATIONS 
Advance Maintenance Dredging 


130. The following recommendations on advance maintenance dredging are 
based not only on Plan 1 results from Model B but also on all other available 
information such as field surveys and field experience. It must be emphasized 
that these recommendations are for an average year including normal storms and 
do not allow for the effects of abnormal storms such as hurricanes and tropi- 
cal storms. Since the Navy wants a minimum clear depth of 46 ft mlw always 
and since it plans to dredge the channel only once a year, the deposition 
rates by zones as well as the deposition rates predicted by Model B for indi- 
vidual cells (Figure 29) have been taken into account in making these 


recommendations. 
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Figure 29. Comparison of computed erosion/deposition rates by cells for 
base and Plan 1 
131. For reasons previously mentioned, the high local deposition rates 
predicted by Model B for base and Plan 1 in some reaches of the channel be- 
tween sta 325+00 and sta 399+74 are suspect because up-to-date bathymetric 
data were not available for model calculations and there are no corroborating 
field data for such high rates. On the other hand, field surveys taken in 


April 1984 and December 1984 which covered the channel between these stations 


Ud 


and which became available after Model B was run, seem to indicate erosion 
rates of the order of 0.8 ft/year or less between sta 325+57 and sta 361+74 
and deposition rates of the order of 1.4 ft/year or less between sta 361+74 
and sta 399+74. In view of the uncertainty on sedimentation rates in the 
length of channel between sta 325+00 and sta 399+74, and because the existing 
depths in this reach are generally higher than 49 ft mlw, an advance mainte- 
mance depth of 3 ft is recommended in this reach. 

132. For convenience in dredging, the channel was divided into reaches 
of at least 2,000-ft lengths at the suggestion of CESAJ. Table 5 lists the 
various reaches of Plan 1 entrance and offshore channels where shoaling is 
expected, estimates of deposition rates (rounded to 0.1 ft/year), and recommen- 
dations for advance maintenance depths (rounded generally to whole feet). If 
the length of channel between sta -97+76 and sta 325+00 is considered and only 
the rectangular portion of the planned channel cross section is taken, the 
total dredging volume for advance maintenance in accordance with the recommen- 
dations shown in Table 5 represents a savings of approximately 630,000 cu yd, 
or nearly 27 percent, compared to the dredging volume for a channel with 3-ft 
advance maintenance throughout this reach. 

133. The recommendations given in Table 5 do not take into account the 
long-term economic advantages of providing greater advance maintenance depths 
and dredging less frequently than once a year, especially in the offshore 
areas, in view of the high cost of mobilization of dredging plant. This issue 
should be explored before a final decision is made on advance maintenance 
depths. 

134. From the geologic sections provided by CESAJ, rock seems to be 
present at depths of 40 to 54 ft mlw between sta 234+00 and 260+00. Two of 
the reaches where large advance maintenance depths of the order of 7 to 9 ft 
have been recommended are in this general area. This is the area just outside 
of the jetty tips and just interior to the jetties. Severe deposition prob- 
lems have been experienced in this general area at present because of material 
moving from the shoals on either side of the channel into the channel. Gener- 
ally, the highest deposition rates have been observed in the northernmost 
quadrant of the channel. In view of the difficulty and expense of dredging in 
rock and problems that may be experienced with large overdepth dredging, it is 
suggested that overwidth dredging be explored as an alternative to overdepth 


dredging in this area. For instance, in addition to providing a reasonable 
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advance maintenance depth, the channel may be widened by 100 to 125 ft (total 
width of channel equals 600 to 625 ft) in this reach. Overwidth dredging may 
be considered also as an alternative in other reaches where rock may be 


present. 


Future Testing 


135. Model B results for Plan 1 provided in this report have been ob- 
tained by testing a 500-ft-wide channel with a project depth of 46 ft mlw and 
advance maintenance depth of 3 ft throughout. Once the channel design is 
finalized, it is recommended that the final design be tested in Model B with 
the latest available bathymetry in and around the channel so that maintenance 
dredging requirements can be determined more accurately corresponding to the 
final channel design. 

136. Model B has been used in this study to estimate average yearly 
erosion/deposition rates in the entrance and offshore channels under average 
wave conditions, excluding abnormal storms. These estimates are good for pre- 
dicting long-term average maintenance dredging requirements. However, the 
effect of severe storms such as hurricanes and tropical storms on shoaling can 
be quite dramatic. So it is recommended that shoaling of the planned naviga- 
tion channel under severe storm conditions be investigated since estimates of 
shoaling volumes can be used in channel design as well as in advance planning 
for emergency mobilization of the necessary dredging plant to keep the channel 
open. This task can be performed using Model B and the storm surge modeling 


capability of WIFM. 
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PART VII: SUMMARY AND CONCLUSIONS 


137. To study the effects of proposed modification of the exterior 
channels of St. Marys Inlet (the ocean entrance to Kings Bay Naval Submarine 
Base) on coastal processes, the CIP system of numerical models of CEWES was 
employed. The system included models for tides, waves, wave-induced currents, 
and sediment transport. The system together with two computational grids 
developed for the study was called Model B. 

138. Model B was used to study existing (base) conditions as well as 
planned conditions. Plan 1 is to (a) widen the navigation channel by 100 ft 
on the north side so the total width becomes 500 ft, (b) deepen the channel to 
-49 ft mlw (46-ft project depth plus 3-ft advance maintenance) with side 
slopes of 3H:1V, and (c) extend the channel on the ocean side with a 20-deg 
bend to the south at sta -97+76. It is assumed also that the landward 
1,000 ft of the south jetty is made sand-tight for Plan l. 

139. The tidal model was verified using the field data of 10 November 
1982. This was achieved by forcing the model with measured tidal elevations 
and matching observed velocities at ranges in the inlet, Cumberland Sound, and 
St. Marys River. There was good agreement. 

140. The average year's wave climate for the study area was obtained 
from WESWIS, on the basis of 20-year hindcast data. The data set included 
normal storms but not hurricanes and tropical storms. This was used in run- 
ning the wave and wave-induced current models. 

141. The sediment transport model determined noncohesive sediment 
(sand) transport in the study area, under the combined action of tides, waves 
and wave-induced currents. It considered a mean tide and the average year's 
wave climate. 

142. The sediment transport model was verified by comparing computed 
erosion/deposition rates in the navigation channel with those obtained from 
field surveys taken by CESAJ during 1980-81. There was good agreement with 
respect to both trends and magnitudes. 

143. While all four models were run for base conditions, only the tide 
model and the sediment transport model were run for plan conditions (Plan 1) 
to meet the urgent need for model results. (Plan conditions were expected to 
influence the tide and sediment transport much more than the waves and 


wave-induced currents. ) 


80 


144. The effects of Plan 1 on tidal currents were mainly local and 
caused by sand-tightening of the south jetty. Velocities at the end of the 
jetties and at tide Gage 1 increased by approximately 10 percent. There were 
no significant changes in velocities at the throat of the inlet, including 
Range 1 and the Fort Clinch area. 

145. Model B predicts an increase in deposition and erosion rates be- 
tween sta -97+76 and sta 325+00 from base to Plan 1. For the reach of channel 
between sta -—80+00 and sta 325+00, the predicted yearly shoaling volumes are 
475,000 and 788,000 cu yd/year for base and Plan 1, respectively, or an in- 
crease of 66 percent for Plan 1. 

146. On the basis of Model B results and all other available informa- 
tion, recommendations on advance maintenance dredging were made for different 
reaches of the navigation channel (Table 5). 

147. For the length of channel between sta -97+76 and sta 325+00, if 
only the rectangular portion of the planned channel cross section is consid- 
ered, the total dredging volume for advance maintenance in accordance with 
Model B recommendations represents a savings of approximately 630,000 cu yd or 
nearly 27 percent compared to the dredging volume for a channel with 3-ft 
advance maintenance throughout according to Plan 1. 

148. In summary, the study successfully accomplished all the study 
objectives, as set forth in paragraph 3, "Purpose," except for the determina- 
tion of waves and wave-induced currents under plan conditions. The numerical 
models for these processes were not rerun, as originally planned, in order to 


meet the urgent needs of the sponsor. 
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Table 1 
Numerical Gages Used in WIFM 


Numerd caliGageinOe spe ay ete ti men Ly MCCA Sem aries een aan 
1 Prototype tide Gage 1 (south spit) 
5 Endeco velocity Gage 2 (main channel) 
10 Range survey Gage 1-A 
iil Range survey Gage 1-B 
2, Range survey Gage 1-C 
25 Ocean end of jetties—channel 
Bi Fort Clinch 
28 South channel (Amelia Island) 


29 Channel to Cumberland 


Table 2 
Sample of WESWIS Data for St. Marys Inlet 
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Table 3 


Incident Wave Conditions for St. Marys Inlet 


DERE EIEN DEES) HEIGST (FT) PERIOD (SECS) FESCENT CUM. FRED. OCCUR. 
2. ei 2 043 


1 -? 22 3.20 4.31 0434 
2 -72.50 2.26 5.39 1,917 06228 
3 -72,50 4.40 6.30 235 18463 
4 -62.5 82 2.50 2.102 OSE65 
5 -52.50 2.46 4.5 £1233 “09793 
i -62.50 4.10 6.20 58 30378 
? -62.50 5.74 7.09 220 30538 
8 -32.50 132 2.40 $1853 32864 
9 -22.50 2.46 4°70 11315 13776 
10 =59159 4:10 5:70 404 141g 
it -=2,50 5.74 5.80 132 gai2 
12 -£2.50 7.38 7.40 1032 14404 
13 -43/59 182 5120 2.733 117202 
14 ~£2,50 2.46 5.50 1.635 19897 
15 -42,50 4.10 5.79 444 19344 
16 -£2,59 5.74 6.20 116 119457 
17 -62,50 7.33 2.20 1 13572 
18 -32.59 22 7.20 9.409 263981 
19 -22, 2.48 5.30 11235 30376 
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22 -32,59 7.38 7.20 029 131077 
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3 -13150 3182 7.40 $99 41305 
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53 27.50 3.74 6.20 370 178366 
53 27.59 7.38 6.90 25 172220 
60 27.30 9:02 7.40 208 72428 
61 27.50 10.66 7.50 155 172593 
62 37:50 182 2.69 73321 
63 37.50 2.46 4.20 11244 174565 
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63 37.50 5.74 6.20 1498 175969 
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77 57.50 5.74 6.10 070 184948 
23 67.59 22 4.90 4.051 22399 
? 67.50 2.46 5.79 ‘624 123693 
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Table 4 
Details of CESAJ Examination Surveys 


Sue a nLiclustvclDatessaa ln lll ml NiSUrveyedlistattons 
1 21-25 Nov 80 130+00 to 325+00 
2 8-9 Dec 80 -80+00 to 130+00 
3 31 Mar-13 Apr 81 130+00 to 325+00 
4 13 Mar-13 Apr 81 -80+00 to 130+00 
5 8 Jun 81 130+00 to 325+00 
6 8-10 Jun 81 -80+00 to 130+00 
7 14-18 Dec 81 -80+00 to 325+00 
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APPENDIX A: NOTATION 


Wave amplitude function 


Ratio of volume of solids to total volume of sediment 
=l-p 


Mapping constants for region p in x-direction 
Mapping constants for region q in y-direction 
Area of cell 

Drag coefficient, wave celerity 

Coefficient 

Wave group velocity 

Chezy coefficient 

Local total water depth 


Sediment diameter exceeded in size by 65 percent (by weight) 
of sediment sample, energy dissipation term 


Wave energy density = peH-/8 

Coriolis parameter, drag coefficient 

Terms representing external forces 

Wave friction factor with D as bed roughness 
Acceleration due to gravity 

Local still-water depth 

Wave height 

Wave height in deep water 


Unit vectors in the x- and y-directions 

Local immersed weight longshore transport rate 
Total immersed weight longshore transport rate 
Number of water cells within the surf zone 
Wave number 

Empirical coefficient 

Empirical coefficient 

Bottom slope 

Maximum values of cell indices for RCPWAVE 


Ratio of group velocity to wave celerity = c /c , Manning's 
roughness o 


Porosity of sediment 


Local volumetric sediment transport rate due to currents 


Al 


Local volumetric sediment transport rate 


Local volumetric sediment transport rates in x- and 
y-directions 


Rate of water volume change due to rainfall or evaporation 


Arbitrary variable, mass density of sediment relative to that 
of fluid (specific gravity of sediment), wave phase function 


Radiation stresses 

Time 

Wave period 

Wave orbital velocity at the bottom 


Time average of the absolute value of the wave orbital 
velocity at bottom 


Tidal velocity components 

Velocity components due to wave-induced currents 
Total velocity component = u + U 

Shear (friction) velocity 

Longshore velocity 

Total velocity component =v +V 
Coordinates in real space 

Width of surf zone 

Dimensionless grain diameter 
Coordinates in computational space 
Breaking index = H/h 

Proportionality coefficient 

Centered difference operators 
Time-step 

Cell dimensions in real space 

Cell dimensions in computational space 
Eddy viscosity for tidal model 

Eddy viscosities in x- and y-directions 
Bed elevation 

Tidal elevation above datum 

Mean free surface displacement (setup) 


Hydrostatic water elevation due to atmospheric pressure 
differences 


Angle of wave propagation 


Wave direction in deep water 


A2 


6 Contour angle 


: Rate of energy dissipation coefficient 

kK Refraction coefficient 

Se Shoaling coefficient 

Ue» n Grid expansion coefficients 

v Kinematic viscosity of fluid 

T 3.14159... 

9) Mass density of sea water 

D . Mass density of solids 

fo) Wave angular frequency = 21/T 

cA Bed shear stress 

eee bes, Bottom friction stresses in x- and y-directions 

toed Lateral shear stress due to turbulent mixing 

> Complex velocity potential for wave 

Superscripts 

k-1 Previous time level 

k Present time level 
k+1 Next time level 

x Intermediate time level 

Subscripts 

b At breaking 

s Stable level of a variable 

ft Partial derivative with respect to time 
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mlw 
msl 
mwl 
NGVD 
NOS 
OICC 
RCPWAVE 
SC 

swl 
USGS 
WESWIS 
WIFM 


APPENDIX B: ABBREVIATIONS AND ACRONYMS 


alternating-direction-implicit 

Coastal Engineering Research Center 

US Army Engineer District, Jacksonville 

US Army Engineer Waterways Experiment Station 
Coastal and Inlet Processes 

Hydraulics Laboratory 

mean low water 

mean sea level 

mean water level 

National Geodetic Vertical Datum 

National Ocean Service 

Officer In Charge of Construction 

Regional Coastal Processes Wave Propagation Model 
stabilizing-correction 

still-water level 

United States Geological Survey 

WES Wave Information Study 

WES Implicit Flooding Model 
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